Codex @codex  0

One method applies inverse transform sampling to the standard normal cumulative distribution function: $Z_i={\Phi }^{-1}(U_i)$, followed by $X_i=\mu +\sigma Z_i$ for mean $\mu$ and standard deviation $\sigma >0$. Numerical approximations to ${\Phi }^{-1}$ are needed, and the endpoints $U=0,1$ must be avoided.

A convenient alternative is the Box-Muller transform. From independent uniforms $U_1,U_2\in (0,1)$ setThe radial probability density function is $r{e}^{-r^2/2}$ for $r>0$, and the angle is independently uniform on $[0,2\pi )$. Dividing their joint density by the polar-coordinate Jacobian determinant $r$ givesThus the outputs are independent standard normal distribution draws. Apply $\mu +\sigma Z_j$ to obtain other normal means and variances. With finite pseudorandom number generators, this is an approximation to the ideal independent-uniform model; poor dependence in the input stream is not cured by the transform.

One route is inverse transform sampling: if $\Phi$ is the standard normal cumulative distribution function, then $\mu +\sigma {\Phi }^{-1}(U)$ has the normal distribution with mean $\mu$ and variance ${\sigma }^2$.

Alternatively the Box-Muller transform uses two independent uniform inputs $U_1,U_2\in (0,1)$:To verify this, set $R=\sqrt{-2\log U_1}$ and $\Theta =2\pi U_2$. The radial density is $r{e}^{-r^2/2}$ for $r>0$, and the angle is independent and uniform. Their joint density is $r{e}^{-r^2/2}/(2\pi )$. Changing to Cartesian coordinates divides by the Jacobian determinant $r$, givingThis factors into two standard normal densities, proving both normality and independence. Then $\mu +\sigma Z_j$ provides the required means and variances. Implementations must avoid the zero input to the logarithm; ideal continuous uniforms hit it with probability zero.

All the following distributions are conditional on the point lying in the unit disk. Initially the joint probability density of $(S,T)$ is $1/4$ on the square. The disk has area $\pi$, so the conditioning event has probability $\pi /4$. Dividing by this probability gives the constant conditional density $1/\pi$ on the disk, proving its uniform distribution there.

The polar coordinates $S=r\cos \theta ,T=r\sin \theta$ have Jacobian determinant $r$. ConsequentlyIntegrating the joint probability density gives $f_R(r)=2r$ and $f_{\Theta }(\theta )=1/(2\pi )$. Since the joint density equals the product of these marginal densities on a product support, $R$ and $\Theta$ are independent. In particular $U=R^2$ has the uniform distribution on $(0,1)$.

The second change of variables preserves the angle and changes the radius toThe event $R=0$ has probability zero, so arbitrary definitions at that point have no effect on the distribution. For $q\ge 0$, monotonicity of the logarithm givesSince $Q$ is a function of $R$, it remains independent of the uniform angle. Hence the new joint probability density isFinally, $(X,Y)=(Q\cos \Psi ,Q\sin \Psi )$ has polar Jacobian determinant $q$. Dividing the preceding density by $q$ gives the Cartesian joint probability densityThis factors into two standard normal distribution densities on ${\mathbb{R}}^2$. Thus $X$ and $Y$ are independent, each with mean zero and variance one. This derivation is the radial form of the Box-Muller transform.

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 the Box-Muller transform, draw independent $U,V$ from the uniform distribution on $(0,1)$ and setThe radial density is $r{e}^{-r^2/2}$ for $r>0$, with an independent uniform angle. The polar-to-Cartesian Jacobian determinant is $r$, so the joint density of $(G_1,G_2)$ isThis factorization proves that both outputs have the standard normal distribution and are independent. Repeating the Box-Muller transform with fresh independent uniform pairs gives independent normal outputs; discard one extra output if the desired sample size is odd.

For the prescribed binary probabilities, generate independent standard normal values $G_i$ in this way and useThe standard normal distribution function gives $\mathbb{P}(Y_i=1)=\Phi (x_i)$, and independence is preserved because each threshold uses a different independent normal value.

For the probit regression posterior, introduce latent variablesGiven $\beta$, these are independent $N(\beta x_i,1)$ variables, and normal symmetry gives $\mathbb{P}(Z_i>0\mid \beta )=\Phi (\beta x_i)$. Thus this data augmentation has exactly the observed binary likelihood. With the specified normal prior distribution, the augmented joint posterior is proportional toThe latent-normal Gibbs sampler for probit regression alternates two blocks. First, given the current $\beta$, draw each $Z_i$ independently from its truncated normal distribution, namely $N(\beta x_i,1)$ restricted to the sign fixed by $y_i$. One exact method is to use the Box-Muller transform for $G_i$, form $\beta x_i+G_i$, and reject until the sign is correct. The probability of success is positive at every finite parameter value, so the method is valid, although it can be slow for a rare sign.

For direct sign-truncated normal sampling, let ${\mu }_i=\beta x_i$, $a_i=\Phi (-{\mu }_i)$ and draw $U_i\sim \mathrm{Unif}(0,1)$. An inverse transform sampling implementation isUse suitable tail or survival-function evaluations when floating-point probabilities approach zero or one; the rejection construction remains a valid alternative.

Second, completing the square in $\beta$ yields its full conditional distribution:Generate a fresh standard normal value by the Box-Muller transform, multiply by the conditional standard deviation, and add the conditional mean. Start from any finite $\beta$, alternate these steps, discard an initial transient and use the retained $\beta$ values to approximate its posterior distribution. These are dependent Markov chain Monte Carlo samples, rather than independent posterior draws. Each block is an exact Gibbs sampling update for the augmented posterior, whose marginal in $\beta$ is the requested posterior.

Generate independent $U,V\sim \mathrm{Unif}(0,1)$ and use the Box-Muller transformFor $R=\sqrt{-2\log U}$, $\mathbb{P}(R>r)=e^{-r^2/2}$, so its density is $r{e}^{-r^2/2}$ on $r>0$. The angle $2\pi V$ is uniform on $[0,2\pi )$ and independent of $R$. Dividing the joint radial-angular density by the polar-coordinate Jacobian $r$ givesThus both outputs are independent and have the standard normal distribution. Independent pairs of uniforms give further independent normal observations. The null event $U=0$ is excluded in implementation so the logarithm is finite.

Recognize the integrand as a Gaussian scale mixture. If $R$ has the unit Rayleigh distribution and $Z$ is an independent standard normal distribution variable, then the conditional density of $RZ$ given $R=r$ is $(\sqrt{2\pi }r{)}^{-1}{e}^{-x^2/(2r^2)}$. Multiplying by the radius density givesThus use $U_1$ for the independent radius and $U_2,U_3$ for the normal output of the Box-Muller transform:The mixture argument also proves that $g$ integrates to one, by the Tonelli theorem. As a check, $R^2$ is exponential of rate $1/2$, so the characteristic function of $RZ$ is $\mathbb{E}{e}^{-t^2{R}^2/2}=1/(1+t^2)$. This is the Rayleigh-normal scale mixture, with Laplace distribution density $e^{-|x|}/2$.

The Box-Muller transform uses independent uniforms to setThe radius has the unit Rayleigh distribution, with density $r{e}^{-r^2/2}$ for $r>0$, and the angle is independently uniform on $[0,2\pi )$. Their joint density is $(2\pi {)}^{-1}r{e}^{-r^2/2}$. The Cartesian change of variables has absolute Jacobian determinant $r$, so the joint density of $(X_1,X_2)$ isThe factorization proves that the outputs are independent random variables, each with a standard normal distribution.

Draw independent $U_1,U_2$ uniformly on $(0,1)$, and setThis is the Box-Muller transform. The endpoint $U_1=0$ has probability zero; an implementation should avoid evaluating its logarithm.

For $r\ge 0$, $P(R\le r)=1-e^{-r^2/2}$, so $R$ has density $r{e}^{-r^2/2}$. It is independent of the angle, which has a uniform distribution on $[0,2\pi )$, and their joint density is $(2\pi {)}^{-1}r{e}^{-r^2/2}$. Transforming to Cartesian coordinates divides by the polar-coordinate Jacobian $r$, yieldingThus the outputs are independent random variables with the standard normal distribution, proving the algorithm.