# Generation of random numbers from the multivariate normal distribution in spherical coordinates

Can any one point me to a research article or post the actual algorithm to generate random numbers from the multivariate normal distribution in spherical coordinates?

I need both $N_p(0,I_p)$ , the standard multivariate normal case in spherical coordinates and $N_p(0,\Sigma)$ where $\Sigma \neq I_p$

• You need to explain first what you want to do. Are you looking to draw samples that are normal distributed in Euclidean x-y-z space and convert the result into spherical coordinates? Or are you interested in drawing samples from a normal distribution in r-phi-theta space? The two are fundamentally different and, in fact, the latter doesn't make much sense because phi and theta are periodic variables and so having a normal distribution in these variables makes not much sense because $p(\phi)$ must equal $p(\phi+2\pi)$ but this is of course inconsistent with the normal distribution.
– Wolfgang Bangerth
Jul 19, 2012 at 22:50
• @Wolfgang Bangerth , i think daek's method is what i need , but i'm not 100% sure yet math.wsu.edu/faculty/genz/papers/mvncmp/mvncmp.html Jul 22, 2012 at 16:38
• THe problem is poorly posed. Do you need $N_p(0,I_p)$? Do you need $N_p(0,\Sigma)$ with $\Sigma \neq I_p$? Do you need $N_p(\mu,\Sigma)$? Aug 6, 2012 at 18:47
• I do not understand the emphasis on spherical coordinates. Is there anything the matter with generating a variable from $N(0,\Sigma)$ and converting it to spherical coordinates?
– whuber
Aug 7, 2012 at 16:09
• Given that even the input to the $\mathcal{N}_p(\mu,\Sigma)$ problem is $O(p^2)$ and the conversion from Cartesian to spherical coordinates is $O(p)$, the time needed for conversion is inconsequential except for small $p$.
– whuber
Aug 8, 2012 at 15:23

If you need $N_p(0,I_p)$, where $p$ is the dimension, then in spherical coordinates, $r\sim \chi^2_p = \Gamma(p/2,1/2)$, and the angles are all independent of each other and the length, with uniform distribution on their respective ranges $[0,\pi)$ for the first one (the one that goes with the cosine only), $[0,2\pi)$ for the remaining ones. If $p$ is even, $p=2m$, you can generate $r_p = \sum_{k=1}^{m} (-\frac12) \ln U_k$, $U_k \sim \mbox{i.i.d. } U[0,1]$. If $p$ is odd, $p=2m+1$, then you need to add another square of a normal to this (see above how to generate them), $r_p = r_{p-1} + z^2$, $z\sim N(0,1)$.
• can you show the proof for $r\sim \chi^2_p = \Gamma(p/2,1/2)$ ? Aug 11, 2012 at 19:31