Generating samples from a Beta(2,2) distribution

Question

I'm looking for a convenient way to generate $\text{Beta}(2,2)$ random variables, using independent $\text{Uniform}(0,1)$ random variables and elementary functions. I'd prefer to avoid rejection or iterative methods.

The PDF is $6 x (1-x)$. The CDF is $x^2(3-2x)$.

I know of one method that seems a bit messy: $X_i\overset{iid}{\sim}\text{Uniform}(0,1)$ $\implies Y_i\equiv-\log(X_i)\overset{iid}{\sim}\text{Exp}(1)$ $\implies \frac{Y_1+Y_2}{Y_1+Y_2+Y_3+Y_4}\sim\text{Beta}(2,2)$. This is a bit inefficient, requiring four uniform samples per beta sample, and it uses logarithms.

The CDF doesn't have a convenient inverse, so an inverse CDF transform won't work either.

EDIT: The median of three $\text{Uniform}(0,1)$ random variables will also follow a $\text{Beta}(2,2)$.

Can anyone think of an elegant method to generate these?

Answer 1

The CDF does have a manageable inverse. Consider:

$$ x^2(3-2x)=u, \;x\in[0,1] $$

That is, you want to know where the CDF crosses the level $u$. We simply need a mapping that tells us what $x$ corresponds to a given $u\in[0,1]$, with the condition that the $x$ we found is in $[0,1]$. This is the same as finding the roots of $x^2(3-2x)-u$ for some fixed $u\in[0,1]$. That's a cubic in $x$ and has at most three roots. I found them using WolframAlpha:

http://www.wolframalpha.com/input/?i=Solve%5Bx%5E2(3-2x)-u,x%5D

It turns out that only one of those gives the result we are after with $x\in[0,1]$: $$ x = -\frac{1}{4}(1+i \sqrt{3}) \sqrt[3]{2 \sqrt{u^2-u}-2 u+1}-\frac{1-i \sqrt{3}}{4 \sqrt[3]{2 \sqrt{u^2-u}-2 u+1}}+\frac{1}{2} $$

I know that looks like it will be a complex number, but it works out to a real number on the interval of interest:

http://www.wolframalpha.com/input/?i=Plot%5B-1%2F4+(1%2Bi+sqrt(3))+(2+sqrt(u%5E2-u)-2+u%2B1)%5E(1%2F3)-(1-i+sqrt(3))%2F(4+(2+sqrt(u%5E2-u)-2+u%2B1)%5E(1%2F3))%2B1%2F2,%7Bx,0,1%7D%5D

I'm guessing there is some simplification to write it without any imaginary units, but the above gets the job done.

So, generate a $u\sim U(0,1)$, plug it into the inverse CDF given above, and that will be $\sim \text{Beta}(2,2)$.

Answer 2

A more elegant computation can be performed to obtain the required inverse of the $\operatorname{Beta}(2,2)$ CDF: if $U \sim \operatorname{Uniform}(0,1)$, then $$X = \frac{1}{2} + \sin \left( \frac{1}{3} \tan^{-1} \frac{U-1/2}{\sqrt{U(1-U)}} \right) \sim \operatorname{Beta}(2,2).$$ Is this computationally efficient? No. Using Mathematica, the command

Total[Take[#, 2]]/Total[#] & /@ Partition[-Log /@ RandomReal[{0, 1}, 4*10^6], 4]

generates $10^6$ realizations of $X$ using $4 \times 10^6$ realizations of $U$ according to the method described in your post. On my machine, this takes about $0.7$ seconds. The command

1/2 + Sin[ArcTan[(-1/2 + #1)/Sqrt[(1 - #1) #1]]/3] & /@ RandomReal[{0, 1}, 10^6]

by contrast, takes about $6.2$ seconds, nearly an order of magnitude slower. This is because it is very much faster to generate extra uniform random variates than to perform trigonometric functions on a single variate. But perhaps a different implementation would reduce this gap in computing time.

Answer 3

Here's elegant solution to the problem, drawn from the idea that the median of three $\text{Uniform}(0,1)$ random variables follows a $\text{Beta}(2,2)$ distribution. We can use this without generating three uniform random variables directly.

First off, a formula exists for the joint density of two order statistics within a random sample, as well as the marginal density of each individual statistic. I can't think of a good reference right now, so I'll just show the specific formulas.

For a sample of size $3$ from a $\text{Uniform}(0,1)$, the joint density of the median, $X_{(2)}$, and the maximum, $X_{(3)}$, is $F_{X_{(2)},X_{(3)}}(x,y)=6x$ for $y\in[0,1]$ and $x\in[0,y]$. Likewise, the marginal density of the maximum is $F_{X_{(3)}}(y)=3y^2$ for $y\in[0,1]$.

From this, we can calculate the conditional density of the median, given the maximum. We get $F_{X_{(2)}|X_{(3)}}(x;y)=\frac{2x}{y^2}$ for $x\in[0,y]$. The conditional CDF is $\frac{x^2}{y^2}$ for $x\in[0,y]$. The inverse CDF is $y\sqrt{x}$ for $x\in[0,1]$.

Now, we can generate a random variable with the density $3y^2$ by $Y_1\equiv \sqrt[3]{U_1}$, where $U_1\sim\text{Uniform}(0,1)$. Then, we can generate a second random variable following the conditional density $\frac{2x}{y^2}$ by $Y_2\equiv y\sqrt{U_2}$, where $U_2\sim\text{Uniform}(0,1)$.

We see then that $\sqrt[3]{U_1}\cdot\sqrt{U_2}\sim\text{Beta}(2,2)$.