# Generating Beta distributions with Uniform generators

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 ?

• A good place to start looking for answers to questions of this form ("how do I generate a random variable from a named distribution") is to search for encyclopedia entries about the distribution: typically, they will include information about random generation of values. The Wikipedia article on the Beta distribution is characteristic.
– whuber
Commented Aug 29, 2015 at 16:16
• Luc Devroye discusses methods for this in his book; look it up there (and his book is free to download, too!) Commented Aug 30, 2015 at 1:06
• The book @Guesswhoitis means is at Luc's pages here Commented Aug 30, 2015 at 2:01

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

• Nice illustration (+1), but another reason it's impractical may be that the beta quantile function has to be approximated numerically & I'd bet rejection methods are quicker (because rbeta uses one). Commented Sep 1, 2015 at 16:23
• @Scortchi Hi there! Thanks for commenting. The post was intended (as all my posts, really) to kind of walk myself through the concepts. The hope is that may be someone else is also wresting with a basic understanding of the idea. There is no intention to provide a method. In any event, please feel free to improve the post by editing. Cheers! Commented Sep 1, 2015 at 16:28
• Oh yes, I got that - just thought it might be useful to mention a reason why you don't always use the probability integral transform in practice: the quantile function isn't always very easy to evaluate. Commented Sep 1, 2015 at 16:38
• I believe the point of the original question may have been to help the student understand the connection between Beta distributions and order statistics of the uniform distribution.
– whuber
Commented Nov 2, 2015 at 16:21
• @whuber It's easy to agree with you, but that's how I interpreted it; not as a computational issue. As a student myself, working on the answer helped me go over this concept again - why is the "boring" uniform actually so interesting. Commented Nov 2, 2015 at 16:24

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

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.