Past exam of the mathematics course of the University of Cambridge / 2012 / iii / Paper 38 / 3

To generate a chi-squared variable with one degree of freedom by the Marsaglia polar method, repeatedly draw independent $V,W$ uniformly on $(-1,1)$ until $0<S=V^2+W^2<1$. ReturnThe accepted point is uniform on the unit disk; its squared radius is uniform on $(0,1)$ and independent of angle. The radial rescaling gives a standard normal coordinate $Z$, hence $Y\sim {\chi }_1^2$. The same accepted pair also supplies another independent standard normal coordinate, if needed. The proof stages below follow the four printed sampler labels.

Squaring a variable with the standard normal distribution gives the required chi-squared distribution. Its two inverse images at $y>0$ are $\pm \sqrt{y}$, so the change-of-variable formula givesDraw the branch-selection uniform variable independently of the randomness used for the normal draw. Then the pair $(Y,U)$ has the product law: conditioning on $Y$ leaves $U$ uniform. This independence is exactly what justifies interpreting the subsequent threshold as the conditional probability of selecting the smaller root.

Introduce $t=\mu y/(2\lambda )\ge 0$. Solving the quadratic and using its product of roots givesBoth roots are positive, with $x_{-}\le \mu \le x_{+}$; equality occurs only for $y=0$. For computation, the smaller root has the equivalent cancellation-free formThis is stable root evaluation for inverse Gaussian sampling.

Define $h(x)=\lambda (x-\mu {)}^2/({\mu }^{2}x)$. It decreases from infinity to zero on $(0,\mu )$ and increases from zero to infinity on $(\mu ,\infty )$, sinceThe two roots are precisely its inverse branches. This supplies the Jacobians needed for the reciprocal-root inverse Gaussian sampler.

For $0<x<\mu$, only the smaller-root branch can produce a value near $x$. Conditional on $Y=h(x)$, its selection probability is $\mu /(\mu +x)$. The output density on that interval is thereforeUseThe factors cancel to giveThus the selected lower branch has exactly the lower part of the desired inverse Gaussian density, rather than the unweighted density of the smaller root alone.

For $x>\mu$, the corresponding smaller root is ${\mu }^2/x$. The conditional probability of selecting the larger root isIt is the same weight as a function of the actual selected value $x$. HenceCombining the two open intervals proves the output has the Inverse Gaussian distribution $IG(\mu ,\lambda )$. The event $Y=0$ has probability zero, so no atom occurs at $\mu$; the limiting density there is finite. The algorithm always selects one of its positive roots, so these branch densities exhaust its probability distribution.

For the constructed output, its defining quadratic says $h(X)=Y$ identically. Since the preceding branch calculation proves that its law is $IG(\mu ,\lambda )$, any variable with that law has the same pushforward distribution:One can also verify this chi-squared transform of an inverse Gaussian variable directly. The density contributions from the inverse branches areThe final equality follows from $x_{-}{x}_{+}={\mu }^2$. This confirms the transformed law independently of which reciprocal root was selected.