How to generate from this distribution without inverse in R/Python?

Question

I am working with a distribution with the following density: $$f(x) = - \frac{(\alpha+1)^2 x^\alpha \log(\beta x)}{1-(\alpha + 1)\log(\beta)}$$ and CDF $$\mathbb{P} (X \leq x) = \int_0^x - \frac{(\alpha+1)^2 t^\alpha \log(\beta t)}{1-(\alpha + 1)\log(\beta)} \, dt = \frac{x^{\alpha+1}((\alpha+1)(\log(\beta x))-1)}{(\alpha+1)\log(\beta)-1}$$ with $x \in (0,1), \beta \in (0,1)$ and $\alpha >-1.$ How can I generate random samples from this distribution in Python/R? Which books can I use to learn about the simulation of random variables and random numbers? Any help is appreciated.

Answer 1

A quick and dirty solution is to apply the inverse transform sampling in which the quantile function is computed via numerical inversion. Below is an R implementation along with an example.

# (your density)
myf <- function(x, a, b) {
 num = (a+1)**2 * x**a * log(b*x)
 den = 1 - (a+1)*log(b)
 return(-num/den)
}

# (your cdf)
myF <- function(x, a, b) {
  num = x**(a+1)* ((a+1)*log(b*x)-1)
  den = (a+1) * log(b) - 1
  return(num/den)
}

# quantile function computed numerically
myqf <- function(p, a, b) {
  uniroot(function(x) myF(x,a,b) -p, 
          lower = .Machine$double.eps, 
          upper = 1-.Machine$double.eps)$root
}

# get some uniform draws
u <- runif(1e+5)

# apply the inverse transform
x <- sapply(u, myqf, a = 5, b=0.99)

hist(x,probability = TRUE)
plot(function(x) myf(x,a=5, b=0.99), from=0, to=1,
     col=2, lwd=2, add= TRUE)

Surely there are less dirty solutions, i.e. more computationally efficient methods to handle this problem, one of these is suggested by Galen in the comments. Many other methods can be found in these three outstanding books:

• Monte Carlo Statistical Methods by Robert and Casella, Springer;
• Computational Statistics by Givens and Hoeting, Wiley;
• Non-Uniform Random Variate Generation by Devroye.

Answer 2

Your pdf is the same as a transformation of a variable that follows a truncated gamma distribution.

The gamma distribution has shape $k=2$, rate $\alpha+1$ and is truncated at $-\log(\beta)$.

Transformation rule for probability density functions

If we apply the transformation

$$\begin{array}{} Y &=& -\log(\beta X) \\ X &=& e^{-Y}/\beta \end{array}$$

then the pdf for $Y$ is

$$g(y) = f\left({e^{-y}}/{\beta}\right) \cdot \left| \frac{dx}{dy}\right| = \frac{{\alpha^\prime}^2 \beta^{-\alpha^\prime} }{1-\alpha^\prime\log(\beta)} y e^{-\alpha^\prime y}$$

with $\alpha^\prime = \alpha +1$ and the domain of $y$ is $(-\log(\beta),\infty)$

Computational

A computational demonstration in R is:

alpha = 1
beta = 0.5

### sampling from truncated gamma distribution using inverse transform
p = runif(10^5, pgamma(-log(beta), shape = 2, rate = alpha+1), 1)        
y = qgamma(p, shape = 2, rate = alpha+1)

### draw a histogram
hist(exp(-y)/beta, breaks = seq(0,1,1/40), freq =F)

### add lines with your density function 
x = seq(0,1,0.001)
y = - ((alpha+1)^2 * x^alpha * log(beta * x))/(1-(alpha + 1)*log(beta))
lines(x,y)

Alternative transformation

Alternatively, one can also apply an additional shift

$$Z = Y+\log(\beta)$$

or effectively

$$Z = -\log(X)$$

which is a variable in the range $(0,\infty)$ and has a pdf like

$$h(z) = g(z-\log(\beta)) = \frac{{\alpha^\prime}^2 \beta^{-\alpha^\prime} \beta e^{\alpha} }{1-\alpha^\prime\log(\beta)} (z-\log(\beta)) e^{-\alpha^\prime z}$$

which is a mixture of two gamma distributed variables (or more specifically Erlang distributions).

So one should be able to perform the trick of Ben with the log transform of uniform variables (creating exponential/Erlang distributed variables), but without the need for rejections sampling.

The method will go something like

• Draw to uniform variables $U_1$, $U_2$.
• Depending on the value of $U_1$ being above or below some value $u_c$
  • have some sum like $z= a \log U_1 + b \log U_2 + c$
  • or use $z= d \log U_2$.
• Compute $x = e^{-z}$

One will need to work out the values $a,b,c,d,u_c$ by writing out $h(z)$ more precisely as the mixture of $\alpha \exp(-\alpha z)$ and $\alpha^2 z \exp(-\alpha z)$. And one needs to take care that $U_1$ is used both as variable to decide which part of the mixture is used as well as generating the exponential distributed variable.

Answer 3

It is possible to generate this random variable using a rejection-sampling algorithm using uniform random variables. To do this, let $U_1,U_2 \sim \text{IID U}(0,1)$ and define the random variable:

$$Z = -\frac{\log(U_1) + \log(U_2)}{\alpha+1} \sim \text{Gamma}(2, \alpha+1).$$

Using the transformation shown in the excellent answer by Sextus Empiricus, it can be shown that:

$$\begin{align} \mathbb{P}(X \leqslant x) &= \mathbb{P} \bigg( \frac{e^{-Z}}{\beta} \leqslant x \bigg| Z \geqslant |\log(\beta)| \bigg) \\[6pt] &= \mathbb{P}( Z \geqslant |\log(\beta x)| |Z \geqslant |\log(\beta)|). \\[6pt] \end{align}$$

Hence, we can generate the random variable $X$ using a rejection-sampling method that generates the uniform values $(U_1,U_2)$ until we get a corresponding value $Z \geqslant |\log(\beta)|$, and then take the appropriate transformation to get $X$. Below we code this method as the function rdist and we also code the CDF as pdist.$^\dagger$

rdist <- function(n, alpha, beta) {
  
  #Check input n
  if (!is.vector(n))        stop('Error: Input n should be a single numeric value')
  if (!is.numeric(n))       stop('Error: Input n should be a single numeric value')
  if (length(n) != 1)       stop('Error: Input n should be a single numeric value')
  if (as.integer(n) != n)   stop('Error: Input n should be an integer')
  if (min(n) <= 0)          stop('Error: Input n should be a non-negative integer')
  if (n == 0) { return(numeric(0)) }
  
  #Check input alpha
  if (!is.vector(alpha))    stop('Error: Input alpha should be a single numeric value')
  if (!is.numeric(alpha))   stop('Error: Input alpha should be a single numeric value')
  if (length(alpha) != 1)   stop('Error: Input alpha should be a single numeric value')
  if (min(alpha) <= -1)     stop('Error: Input alpha should be greater than minus-one')
  
  #Check input beta
  if (!is.vector(beta))     stop('Error: Input beta should be a single numeric value')
  if (!is.numeric(beta))    stop('Error: Input beta should be a single numeric value')
  if (length(beta) != 1)    stop('Error: Input beta should be a single numeric value')
  if (min(beta) <= 0)       stop('Error: Input beta should be greater than zero')
  if (min(beta) >= 1)       stop('Error: Input beta should be less than one')
  if ((alpha+1)*log(beta) == 1) stop('Error: Inadmissible parameters')
  
  #Generate pseudo-random values
  x <- rep(0, n)
  min.z <- -log(beta)
  for (i in 1:n) {
    z <- 0
    while (z <= min.z) { u <- runif(2); z <- -sum(log(u))/(alpha+1) }
    x[i] <- exp(-z)/beta }
  
  #Give output
  x }

#####################################################################

pdist <- function(x, alpha, beta, lower.tail = TRUE, log.p = FALSE) {
  
  #Check input x
  if (!is.vector(x))        stop('Error: Input x should be a numeric vector')
  if (!is.numeric(x))       stop('Error: Input x should be a numeric vector')
  if (length(x) == 0) { return(numeric(0)) }
  
  #Check input alpha
  if (!is.vector(alpha))    stop('Error: Input alpha should be a single numeric value')
  if (!is.numeric(alpha))   stop('Error: Input alpha should be a single numeric value')
  if (length(alpha) != 1)   stop('Error: Input alpha should be a single numeric value')
  if (min(alpha) <= -1)     stop('Error: Input alpha should be greater than minus-one')
  
  #Check input beta
  if (!is.vector(beta))     stop('Error: Input beta should be a single numeric value')
  if (!is.numeric(beta))    stop('Error: Input beta should be a single numeric value')
  if (length(beta) != 1)    stop('Error: Input beta should be a single numeric value')
  if (min(beta) <= 0)       stop('Error: Input beta should be greater than zero')
  if (min(beta) >= 1)       stop('Error: Input beta should be less than one')
  if ((alpha+1)*log(beta) == 1) stop('Error: Inadmissible parameters') 
  
  #Check inputs lower.tail and log.p
  if (!is.vector(lower.tail))  stop('Error: Input lower.tail should be a single logical value')
  if (!is.logical(lower.tail)) stop('Error: Input lower.tail should be a single logical value')
  if (length(lower.tail) != 1) stop('Error: Input lower.tail should be a single logical value')
  if (!is.vector(log.p))       stop('Error: Input log.p should be a single logical value')
  if (!is.logical(log.p))      stop('Error: Input log.p should be a single logical value')
  if (length(log.p) != 1)      stop('Error: Input log.p should be a single logical value')
  
  #Generate probabilities
  a <- alpha+1
  n <- length(x)
  PROBS <- rep(0, n)
  for (i in 1:n) {
    xx <- x[i]
    if (xx >= 1) { PROBS[i] <- 1 }
    if ((xx > 0)&(xx < 1)) {
      PROBS[i] <- (xx^a)*(a*log(beta*xx)-1)/(a*log(beta)-1) } }
  if (!lower.tail) { PROBS <- 1-PROBS }
  
  #Give output
  if (log.p) { log(PROBS) } else { PROBS } }

We can now use this function to generate a large number of pseudo-random values. Plotting the resulting ECDF for the values gives a close approximation to the true CDF.

#Set distribution parameters
alpha <- 2
beta  <- 0.6

#Generate values from the distribution
set.seed(1)
n <- 10^6
X <- rdist(n, alpha, beta)

#Plot empirical cumulative distribution function (ECDF) of values
plot(ecdf(X), col = 'blue', main = 'ECDF of pseudo-random values', ylab = 'Proportion')

---

$^\dagger$ Here is alternative code for rdist if you are want to avoid rejection sampling by using the pgamma and rgamma functions.

rdist <- function(n, alpha, beta) {
  
  #Check input n
  if (!is.vector(n))        stop('Error: Input n should be a single numeric value')
  if (!is.numeric(n))       stop('Error: Input n should be a single numeric value')
  if (length(n) != 1)       stop('Error: Input n should be a single numeric value')
  if (as.integer(n) != n)   stop('Error: Input n should be an integer')
  if (min(n) <= 0)          stop('Error: Input n should be a non-negative integer')
  if (n == 0) { return(numeric(0)) }
  
  #Check input alpha
  if (!is.vector(alpha))    stop('Error: Input alpha should be a single numeric value')
  if (!is.numeric(alpha))   stop('Error: Input alpha should be a single numeric value')
  if (length(alpha) != 1)   stop('Error: Input alpha should be a single numeric value')
  if (min(alpha) <= -1)     stop('Error: Input alpha should be greater than minus-one')
  
  #Check input beta
  if (!is.vector(beta))     stop('Error: Input beta should be a single numeric value')
  if (!is.numeric(beta))    stop('Error: Input beta should be a single numeric value')
  if (length(beta) != 1)    stop('Error: Input beta should be a single numeric value')
  if (min(beta) <= 0)       stop('Error: Input beta should be greater than zero')
  if (min(beta) >= 1)       stop('Error: Input beta should be less than one')
  if ((alpha+1)*log(beta) == 1) stop('Error: Inadmissible parameters')
  
  #Generate pseudo-random values
  min.p <- pgamma(-log(beta), shape = 2, rate = alpha+1)
  z <- qgamma(runif(n, min.p, 1),  shape = 2, rate = alpha+1)
  x <- exp(-z)/beta
  
  #Give output
  x }

Answer 4

There are already some great answers here. Another way to do this would be to use MCMC a la the Metropolis Hastings Algorithm. See my implementation below

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta

@np.vectorize
def target_distribution(x):
    
    a = 5
    b = .999
    
    top = -(a+1)**2 * x**a * np.log(b*x)
    bottom = 1 - (a+1)*np.log(b)
    
    return top/bottom



# Initialize the Metropolis-Hastings algorithm with an initial value of x = 0 and a number of iterations
x = 0.5
iterations = 100000

# Initialize an array to store the samples
samples = np.zeros(iterations)

# Run the Metropolis-Hastings algorithm
for i in range(iterations):
    # Sample from the proposal distribution
    x_proposal = beta(a=1, b = 1).rvs()
    
    # Calculate the acceptance ratio
    acceptance_ratio = target_distribution(x_proposal) / target_distribution(x)
    
    # Accept or reject the proposal
    if np.random.uniform() < acceptance_ratio:
        x = x_proposal
    
    # Store the sample
    samples[i] = x

# Plot the samples
plt.hist(samples, bins=50, density=True, alpha=0.5)
x_range = np.linspace(0, 1, 1000)
plt.plot(x_range, target_distribution(x_range))
plt.show()

Answer 5

Though you seem to be willing to avoid inversion, this can be done analytically.

After setting $z:=x^{\alpha+1}$, from $$u(z):=\text{cdf}(z)=z(\gamma\log z+1)=\gamma e^{-1/\gamma}(e^{1/\gamma}z)\log(e^{1/\gamma}z)=\gamma e^{-1/\gamma} W^{-1}(e^{1/\gamma}z)$$

we draw

$$e^{1/\gamma}z(u)=e^{1/\gamma}\text{cdf}^{-1}(u)=W\left(\gamma^{-1} e^{1/\gamma}u\right)$$

where $W$ denotes Lambert's function.

Now it suffices to draw $u$ uniformly in $[0,1]$ and pass it to $\text{cdf}^{-1}(u)$.

Answer 6

I think it is possible to simulate this simply using three independent standard uniform random variables, without any approximation or numerical interpolation or rejection sampling or use of the Lambert W function.

Yves Daoust's answer used $Z= X^{\alpha+1}$ and $\gamma= \frac{1}{{(\alpha+1)\log(\beta)-1}}$ to simplify the cdf for $Z$ to $z(\gamma\log z+1)$ on $(0,1]$ and so the density to $1+\gamma +\gamma \log(z)$. You will have $-1 \le \gamma \le 0$.

This is a mixture distribution for $Z$, with probability $1+\gamma$ of being drawn from a uniform distribution on $(0,1]$ and probability $-\gamma$ of being drawn from a distribution with density $-\log(z)$ on the same support. Somebody previously asked about the latter distribution with an answer suggesting it might be a negative log-gamma distribution and a comment suggesting it might be a product uniform distribution as the distribution of the product of two standard uniform random variables. That is the kind of thing we want to make random samples. The R code is short and simple:

rthisdist <- function(n, alpha=0, beta=1){
   # require alpha > -1 and 
   #         beta positive for logarithm and <= 1; 1E-99 will do
   gam <- 1 / ((alpha+1) * log(beta) - 1) # gamma without overloading  
   uniformA <- runif(n)
   uniformB <- runif(n)
   uniformC <- runif(n)
   z <- ifelse(uniformA < -gam, uniformB, 1) * uniformC
   return(z ^ (1 / (alpha+1)))
   }

Testing this against your CDF gives a convincing match (simulations in blue, your CDF in red):

testalpha <- 2   
testbeta  <- 0.6 
plot.ecdf(rthisdist(10^4, testalpha, testbeta), col="blue")
curve((x^(testalpha+1) * ((testalpha+1) * (log(testbeta * x)) - 1)) / 
      ((testalpha+1) * log(testbeta) - 1), 
      from=0, to=1, add=TRUE, col="red")