Generating Beta distributions with Uniform generators

Question

I can generate as many samples from one or more uniform distribution (0,1) as I wish. How can I use this to generate a beta distribution ?

Answer 1

I will use [R] not so much as a practical answer (rbeta would do the trick), but as an attempt at thinking through the probability integral transform. I hope you are familiar with the code so you can follow, or replicate (if this answers your question).

The idea behind the Probability Integral Transform is that since a $cdf$ monotonically increases in value from $0$ to $1$, applying the $cdf$ function to random values form whichever distribution we may be interested in will on aggregate generate as many results say, between $0.1$ and $0.2$ as from $0.8$ to $0.9$. Now, this is exactly what a $pdf$ of a $U(0,1)$. It follows that if we start with values from a random uniform, $U \sim (0,1)$ instead, and we apply the inverse $cdf$ of the distribution we are aiming at, we'll end up with random values of that distribution.

Let's quickly show it with the queen of the distributions... The Normal $N(0,1)$. We generate $10,000$ random values, and plug them into the $erf$ function, plotting the results:

# Random variable from a normal distribution:
x <- rnorm(1e4)
par(mfrow=c(1,2))
hist(x, col='skyblue', main = "Random Normal")

# When transform by obtaining the cdf (x) will give us a Uniform:
y <- pnorm(x)
hist(y, col='skyblue', main = "CDF(X)")

In your case, we are aiming for $X \sim Beta(\alpha, \beta)$. So let's get started at the end and come up with $10,000$ random values from a $U(0,1)$. We also have to select values for the shape parameters of the $Beta$ distribution. We are not constrained there, so we can select for example, $\alpha=0.5$ and $\beta=0.5$. Now we are ready for the inverse, which is simply the qbeta function:

U <- runif(1e4)
alpha <- 0.5
beta <- 0.5
b_rand <- qbeta(U, alpha, beta)
hist(b_rand, col="skyblue", main = "Inverse U")

Compare this to the shape of the $Beta(\alpha,\beta)$ $pdf$:

x <- seq(0, 1, 0.001)
beta_pdf <- dbeta(x, alpha, beta)
plot(beta_pdf, type ="l", col='skyblue', lwd = 3,
     main = "Beta pdf")

Answer 2

If $U_i$ are iid $U(0,1)$ random variables, then $-log(U_i)$ is exponential distribution.

If $X_i$ are iid $Exp(1)$ random variables, then beta distribution can be derived as $$Y=\frac{\sum_{j=1}^{a}X_j}{\sum_{j=1}^{a+b}X_j}\sim Beta(a,b)\quad (a,b\in \mathbb N^*)$$

Nsim=10^4 #number of random variables
a=3;b=5
U=runif(Nsim)
X=-log(U) #transforms of uniforms
E=rexp(Nsim) #exponentials from R
par(mfrow=c(1,2)) #plots
hist(X,freq=F,main="Exp from Uniform")
hist(E,freq=F,main="Exp from R")

Y=c()
for (i in 1:Nsim) {
  rb=runif(a+b)
  Y=c(Y, sum(-log(rb[1:a]))/sum(-log(rb)))
}
B=rbeta(Nsim, shape1 = a, shape2 = b) #beta from R
par(mfrow=c(1,2)) #plots
hist(Y,freq=F,main="Beta from Uniform")
hist(B,freq=F,main="Beta from R")

The caveat of the first simulation is that $a$ and $b$ must be positive integers. There is another method without this restriction:

If $U$ and $V$ are iid $U(0,1)$ random variables, the distribution of $$\frac{U^{1/a}}{U^{1/a}+V^{1/b}}$$ conditional on $U^{1/a}+V^{1/b}\le 1$, is the $Beta(a; b)$ distribution.

Nsim=10^4 #number of random variables
a=0.5;b=0.5

W=c()
for (i in 1:Nsim) {
  U=runif(1)
  V=runif(1)
  while(U^(1/a)+V^(1/b)>1){
    U=runif(1)
    V=runif(1)
  }
  W=c(W, U^(1/a)/(U^(1/a)+V^(1/b)))
}

B=rbeta(Nsim, shape1 = a, shape2 = b) #beta from R
par(mfrow=c(1,2)) #plots
hist(W,freq=F,main="Beta from Uniform")
hist(B,freq=F,main="Beta from R")

Answer 3

You could use the inverse transform sampling method, which is useful to know about because it's a very general method that is not limited solely to the beta distribution.