I was coding a function to generate a 3D Gaussian distribution. However, I obtain unexpected results.
My idea was to generate these numbers $(x,y,z)$ through the following procedure:
- Generate three uniformly distributed numbers $r_0\in[0,1]$, $\theta\in[0,\pi]$, and $\phi\in[0,2\pi]$;
- Use the transformation $r = \sqrt{-\log r_0}$; and
- Use spherical coordinates as $$\begin{align} x &= r \sin\theta \cos\phi \\ y &= r \sin\theta \sin\phi \\ z &= r \cos\theta \end{align}$$
After running this algorithm several times, the distribution of $z$ is effectively gaussian, but the distributions of $x$ and $y$ are not (see following figure). What am I missing?
By the way, I coded this algorithm with Fortran:
subroutine rand_maxwellian(xx, yy, zz)
implicit none
real, intent(out) :: xx, yy, zz
real :: pi = acos(-1), rr, th, ph
call random_seed()
call random_number(th) ! theta = [0,1]
call random_number(ph) ! phi = [0,1]
call random_number(rr) ! radius = [0,1]
th = pi*th ! theta = [0,pi]
ph = 2.0*pi*ph ! phi = [0,2pi]
rr = sqrt( -log(rr) ) ! radius in a gaussian distribution
xx = rr * sin(th) * cos(ph)
yy = rr * sin(th) * sin(ph)
zz = rr * cos(th)
end subroutine rand_maxwellian
The procedure you use is mimicking the standard 2D way to construct such Gaussians. And indeed $z$ is Gaussian and $\sqrt{x^2+y^2}$ is the abs value of a Gaussian (you may check this numerically). But multiplying the latter with $\cos \phi$ and $\sin \phi$ does not produce Gaussians.