= Solution
For the <Box-Muller transform>, draw <independent> $U,V$ from the <uniform distribution> on $(0,1)$ and set
$$
R=\sqrt{-2\log U},\qquad \Theta=2\pi V,\qquad
\boxed{G_1=R\cos\Theta,\quad G_2=R\sin\Theta.}
$$
The 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)$ is
$$
\frac1{2\pi}\exp\left[-\frac{g_1^2+g_2^2}{2}\right]
=\frac{e^{-g_1^2/2}}{\sqrt{2\pi}}\frac{e^{-g_2^2/2}}{\sqrt{2\pi}}.
$$
This 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 use
$$
\boxed{Y_i=\mathbf1_{\{G_i\leq x_i\}}.}
$$
The <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 variables
$$
Z_i=\beta x_i+G_i,\qquad Y_i=\mathbf1_{\{Z_i>0\}}.
$$
Given $\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 to
$$
\exp\left[-\frac{\beta^2}{2}-\frac12\sum_i(z_i-\beta x_i)^2\right]
\prod_i\mathbf1_{\{z_i>0\text{ if }y_i=1;\ z_i\leq0\text{ if }y_i=0\}}.
$$
The <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\operatorname{Unif}(0,1)$. An <inverse transform sampling> implementation is
$$
\boxed{Z_i=\mu_i+\Phi^{-1}(a_i+(1-a_i)U_i)\quad(y_i=1),\qquad
Z_i=\mu_i+\Phi^{-1}(a_iU_i)\quad(y_i=0).}
$$
Use 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>:
$$
\boxed{\beta\mid z,y\sim N\left(\frac{\sum_i x_i z_i}{1+\sum_i x_i^2},\ \frac1{1+\sum_i x_i^2}\right).}
$$
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.
Back to article page
