Past exam of the mathematics course of the University of Cambridge / 2008 / iii / Paper 47 / 3

A good pseudorandom number generator for Monte Carlo integration should have approximately uniform outputs, a period much longer than the intended run, and negligible detectable dependence both serially and in higher-dimensional tuples. It should be fast, reproducible from a seed, portable enough to reproduce a calculation, and numerically precise enough for the required tail probabilities. Values strictly inside $(0,1)$ avoid singularities in logarithmic inverse transform sampling. Separate simulations also need well-separated streams or seeds, rather than accidentally identical streams. A pseudorandom number generator is deterministic; these properties mean that it imitates independent uniform distributions sufficiently well for the application, not that its outputs are literally independent.

For rate $\lambda >0$, the exponential distribution has cumulative distribution function $F(x)=1-e^{-\lambda x}$ on $x\ge 0$. Inverse transform sampling gives $-\log (1-U)/\lambda$; since $1-U$ is also uniform, an equivalent answer isFor $x\ge 0$, $P(E>x)=P(U<e^{-\lambda x})=e^{-\lambda x}$, verifying the exponential distribution directly.

The centered Laplace distribution with density $(\lambda /2)e^{-\lambda |x|}$ can be obtained by multiplying $E$ by an independent sign taking $1$ and $-1$ with equal probabilities. Alternatively, its inverse transform sampling formula uses just one uniform variate:The first branch inverts the negative half of the cumulative distribution function and the second branch inverts the positive half. Add a location parameter if a noncentered Laplace distribution is wanted.

Correct uniformity is essential for each marginal distribution, and independence is essential for independent samples and for the sign construction. For example, setting every $U_i$ equal to a single uniform $U$ gives uniform marginals but identical transformed draws. Likewise, choosing a sign from the same uniform used for the magnitude can change the Laplace distribution unless the one-uniform inverse above is used.

The first algorithm uses independent normal distributions: generate $Z_1,\dots ,Z_{\nu }$ with law $N(0,1)$ and returnFor example, the Box-Muller transform generates two such normal variates from independent uniforms:To verify this, $R^2=-2\log U$ has density $e^{-s/2}/2$ and $2\pi V$ is an independent uniform angle. Converting the resulting polar density to Cartesian coordinates gives $(2\pi {)}^{-1}{e}^{-(z_1^2+z_2^2)/2}$, which factors into two independent normal distributions. Independent pairs therefore give the required sum, discarding a spare coordinate when $\nu$ is odd. Its moment-generating function is $(1-2t{)}^{-\nu /2}$, that of the chi-squared distribution.

A genuinely different algorithm uses exponential-envelope rejection sampling for a chi-squared variable. Put $a=\nu /2\ge 1$. The target and proposal probability density functions areFor $a>1$, differentiating $\log (f/g)$ gives $(a-1)(1/y-1/\nu )$, so the maximum is at $y=\nu$. The envelope constant and normalized density ratio areThey remain valid at $a=1$, where $M=1$ and $f=g$. Repeatedly draw fresh independent $U,V$, put $Y=-\nu \log U$, and accept whenThe right side is nonpositive by $\log u\le u-1$, so this is a legitimate acceptance probability. The joint probability of proposing and accepting in $dy$ is $g(y)f(y)/(Mg(y))dy=f(y)dy/M$; conditional on acceptance the density is exactly $f$. This proves the rejection sampling algorithm, rather than only identifying the target law.

For $\nu =2$ I prefer the second method: it reduces to a single exponential draw with no rejection. The first method is simple and attractive when normal draws are already available or $\nu$ is small. For larger $\nu$, the second method avoids generating $\nu$ normals: its acceptance probability is $1/M$, asymptotic to $\sqrt{2\pi }/(e\sqrt{a})$ by Stirling's formula. Its expected number of proposals grows like $\sqrt{\nu }$, rather than the linear number of normal draws in the first method. Both are exact ideal-uniform algorithms; actual runtime also depends on the implementation.

For Beta sampling by gamma ratios, use either preceding chi-squared distribution algorithm to generate independent $S\sim {\chi }_{2\alpha }^2$ and $T\sim {\chi }_{2\beta }^2$, and returnWrite $W=S+T$, so $S=WR$ and $T=W(1-R)$. The Jacobian determinant is $W$. The joint probability density function after this change of variables factors asIntegrating out $w$ proves the Beta distribution for $R$.

A new method is Beta sampling by uniform order statistics. Generate $m=\alpha +\beta -1$ independent uniforms, sort them, and return the $\alpha$th order statistic. There are $\alpha -1$ uniforms below the selected value and $\beta -1$ above it. Selecting their labels and placing the selected observation in $dr$ gives densitywhich is the Beta distribution density. The integrality of the shape parameters is what makes this order-statistic construction possible.