12
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ You could try the no u turn sampler if you confirm the likelihood is almost everywhere smooth. $\endgroup$
    – Galen
    Apr 9, 2023 at 15:20
  • $\begingroup$ For completeness, it might be worth clarifying what happens in the case where $\log(\beta) = \tfrac{1}{\alpha+1}$ (where you have denominators of zero). $\endgroup$
    – Ben
    Apr 10, 2023 at 22:13
  • 7
    $\begingroup$ Voting to reopen: I disagree strongly that this question is anything other than "a statistical topic disguised as a coding question." Also this may be a good resource Devroye, L. (1986). Non-Uniform Random Variate Generation. Springer Verlag. $\endgroup$
    – Alexis
    Apr 10, 2023 at 23:09
  • 3
    $\begingroup$ The parameter $\alpha$ can be "swallowed" by the change of variable $z=x^{\alpha+1}$, reducing to $\mathbb P(Z\le z)=z(\gamma\log z+1)$. $\endgroup$ Apr 11, 2023 at 6:47

6 Answers 6

17
$\begingroup$

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)

enter image description here

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:

$\endgroup$
3
  • 1
    $\begingroup$ @SextusEmpiricus, I might be wrong but I read it "the quantile function is not available" which is why I proposed to approximate it numerically. $\endgroup$
    – utobi
    Apr 9, 2023 at 20:21
  • 2
    $\begingroup$ Actually, I think that I can use this method. $\endgroup$ Apr 9, 2023 at 20:25
  • $\begingroup$ @utobi you are right, I read it now as well as 'this distribution without inverse'. I read it before as a typo 'without using inverse transform sampling'. $\endgroup$ Apr 9, 2023 at 21:38
14
$\begingroup$

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)

example of comparison between histogram of simulation and curve of density

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.

$\endgroup$
2
  • 1
    $\begingroup$ This Wikipedia page has more details about the transformation of the pdf when a variable is transformed. $\endgroup$ Apr 9, 2023 at 21:59
  • 1
    $\begingroup$ Great answer (+1). Your alternative transformation is not entirely clear to me, and I'm not convinced it would work without a rejection step, but perhaps seeing is believing. $\endgroup$
    – Ben
    Apr 11, 2023 at 22:32
11
$\begingroup$

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')

enter image description here


$^\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 }
$\endgroup$
7
$\begingroup$

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()

enter image description here

$\endgroup$
1
  • $\begingroup$ Nice! My first thought was to use some Monte Carlo method, but I'm new to this area and still learning. I will read about this method, thank you! $\endgroup$ Apr 10, 2023 at 14:01
7
$\begingroup$

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)$.

$\endgroup$
15
  • $\begingroup$ Nice, +1. Yet, W has to be computed analytically en.m.wikipedia.org/wiki/Lambert_W_function $\endgroup$
    – utobi
    Apr 11, 2023 at 7:47
  • 1
    $\begingroup$ @utobi: you probably mean numerically ? $\endgroup$ Apr 11, 2023 at 7:53
  • $\begingroup$ @utobi The rnorm function in R, and probably derived from something else, also uses a numerical approximation with the use of a high order polynomial. The advantage is that you can use the same polynomial independent of the scale and location. This seems to be the case here as well as the numerical approximation boils down to a single function $W$. Because you only need a single function, independent of the parameters, making a numerical approach should be achievable very efficiently. $\endgroup$ Apr 11, 2023 at 7:57
  • $\begingroup$ @SextusEmpiricus: also, accuracy not being critical for RNG, a coarse precomputed tabulation can be used, I guess. $\endgroup$ Apr 11, 2023 at 8:02
  • 1
    $\begingroup$ God dangint, I LOVE the Lambert W Function, good on you +1 $\endgroup$ Apr 11, 2023 at 18:23
2
$\begingroup$

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")

enter image description here

$\endgroup$
2
  • $\begingroup$ It can be done with two uniform variables. After using the uniformA for the if/else condition, you can reuse it as the uniform variable in the further computation. $\endgroup$ Apr 12, 2023 at 14:53
  • $\begingroup$ @SextusEmpiricus - yes indeed, though I think you would have to divide uniformA by $-\gamma$ to replace uniformB and that might make it slightly harder to see what was going on. Another simplification could be to change the sign of $\gamma$ so it was $\frac{1}{1-(\alpha+1)\log(\beta)} \in [0,1]$ and directly represented a probability in the mixture distribution. $\endgroup$
    – Henry
    Apr 12, 2023 at 14:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.