# Monte Carlo simulation on sphere: unbiased random steps

Im doing a Metropolis Monte Carlo simulation with particles on a sphere and have a question concerning the random movement in a given time step.

I understand that to obtain a uniform distribution of random points on the sphere to begin with it is not enough to use the naive simplest way (use spherical coordinates with constant R and pick random angles theta and phi), but one has to for example use one of these methods: http://mathworld.wolfram.com/SpherePointPicking.html

Looking at another code for a Monte Carlo on a sphere I see a fairly complicated way to generate random moves: pick a random particle, calculate the rotation matrix moving it to the north pole, find a random cartesian vector less than a certain length and move it to the north pole, normalize the cartesian vector and then rotate it back to the vicinity of the original particle position.

This is all to get an unbiased new random position. I dont understand the rationale completely although I suspect it is connected to detailed balance. But my feeling is that there should be an easier (i.e. faster) way to do this. Actually, intuitively I feel that in this case it is ok to find two small random angles theta and phi and add them to the position of the particle - or would this be a mistake?

• Please, dont' cross post stackoverflow.com/questions/12900861/… – leonbloy Oct 15 '12 at 19:20
• Sorry, didnt know that was frowned upon. Actually someone in the programming forum suggested I try posting the question here. – jorgen Oct 15 '12 at 20:18

## 2 Answers

$\newcommand{\+}{^{\dagger}}% \newcommand{\angles}[1]{\left\langle #1 \right\rangle}% \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}% \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}% \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}% \newcommand{\dd}{{\rm d}}% \newcommand{\ds}[1]{\displaystyle{#1}}% \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}% \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}% \newcommand{\fermi}{\,{\rm f}}% \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}% \newcommand{\half}{{1 \over 2}}% \newcommand{\ic}{{\rm i}}% \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow}% \newcommand{\isdiv}{\,\left.\right\vert\,}% \newcommand{\ket}[1]{\left\vert #1\right\rangle}% \newcommand{\ol}[1]{\overline{#1}}% \newcommand{\pars}[1]{\left( #1 \right)}% \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}}% \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}% \newcommand{\sech}{\,{\rm sech}}% \newcommand{\sgn}{\,{\rm sgn}}% \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}}% \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ Let's assume we have a mechanism to generate random numbers in $\left[0,1\right)$ and ${\rm P}\pars{\Omega_{\vec{r}}}$ is a distribution function for random points in a sphere of radius $a > 0$. $\Omega_{\vec{r}}$ is the solid angle. In this case, ${\rm P}\pars{\Omega_{\vec{r}}}$ is, indeed, $\Omega_{\vec{r}}\,\,$-independent: $$1 = \int{\rm P}\pars{\Omega_{\vec{r}}}\,\dd\Omega_{\vec{r}} = {\rm P}\pars{\Omega_{\vec{r}}}\int\dd\Omega_{\vec{r}} = {\rm P}\pars{\vec{r}}\pars{4\pi} \quad\imp\quad{\rm P}\pars{\vec{r}} = {1 \over 4\pi}$$ Then, $$1 = \int_{0}^{\pi}\half\,\sin\pars{\theta}\,\dd\theta\int_{0}^{2\pi} \,{\dd\phi \over 2\pi}$$ We can generate random numbers $\xi_{\theta}$ and $\xi_{\phi}$ such that: $$\bracks{~\half\,\sin\pars{\theta}\,\dd\theta = \dd\xi_{\theta}\,, \quad\xi_{0} = 0 \imp \theta = 0~}\ \mbox{and}\ \bracks{~{\dd\phi \over 2\pi} = \dd\xi_{\phi}\,,\quad\xi_{0} = 0 \imp \phi = 0~}$$ Those relations yield: $\ds{\half\bracks{-\cos\pars{\theta} + 1} = \xi_{\theta}}$ $\ds{\pars{~\mbox{or/and}\ \sin\pars{\theta/2} = \root{\xi_{\theta}}~}}$ and $\ds{\phi = 2\pi\,\xi_{\phi}}$: $$\color{#0000ff}{\large% \theta = 2\arcsin\pars{\root{\xi_{\theta}}}\,,\qquad \phi = 2\pi\xi_{\phi}}$$ As we mentioned above, $\xi_{\theta}$ and $\xi_{\phi}$ are uniformly distributed in $\left[0, 1\right)$.

For a sphere of radio $a$ the random points are given by: $$\left\lbrace% \begin{array}{rcl} \color{#0000ff}{\large x} & = & a\sin\pars{\theta}\cos\pars{\phi} = \color{#0000ff}{\large 2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}\cos\pars{2\pi\xi_{\phi}}} \\ \color{#0000ff}{\large y} & = & a\sin\pars{\theta}\sin\pars{\phi} = \color{#0000ff}{\large 2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}\sin\pars{2\pi\xi_{\phi}}} \\ \color{#0000ff}{\large z} & = & a\cos\pars{\theta} = \color{#0000ff}{\large a\pars{1 - 2\xi_{\theta}}} \end{array}\right.$$

I don't clearly see what the second paragraph (uniform distribution) has to do with the rest (unbiased random steps). Anyway: in classical Montecarlo the choosing of displacement is not critical, typically one is only concerned with:

1. sampling the full configuration space (this typicaly means that the displacement is random enough in "direction", so that any configuration is a priori reachable)

2. having a "nice" probability (not almost one, nor almost zero) of move acceptance (this typically means that the average displacement is not very large, nor very small)

In your case, I'd try the simplest alternative: for a selected point choose a (random or fixed) displacement, and choose a random direction (a uniform angle). The "fairly complicated" procedure you describe seem to be just a implementation of this approach.

• Thanks for the reply! According to the code I described, the method is intended to "move the point in an unbiased way" - I was assuming this had something to do with detailed balance (the probability of being in state 1 times the probability to move from 1 to 2 equals the opposite). I guess what I was trying to ask was whether I run into trouble if I, instead of doing this procedure, simply choose a phi and theta of random (different) lengths between zero and some maximum, where theta and phi are the angles of the standard spherical coordinate system, and increment the angles by this? – jorgen Oct 15 '12 at 18:51
• The point of that code is to jump to any point in a neighborhood with equal probability. That's a not strictly necessary but nice-to-have property (and you don't get that by just summing two angles to the plain polar coordinates of your point). – leonbloy Oct 15 '12 at 18:57
• Thanks! That seems like a lot of computer time used for something that is not necessary.. Perhaps convergence etc. is a lot faster this way. – jorgen Oct 15 '12 at 19:16
• @jorgen your approach (make small deltas to phi and theta) will lead to vastly different steps near the pole than near the equator, and that could very easily actually bias your results. So yes, you could very easily run into trouble. – Steven Stadnicki Sep 12 '13 at 6:31
• @jorgen Worse, if you do run into trouble - that is, if your simulator gives inaccurate results - then you might not notice the bias that you've introduced. For purposes of simulation like what you're doing it's probably moot, but it's a good idea to try and do the Right Thing when you can. Furthermore, it's not nearly so complicated as you're making it out to be; in fact, it might well be cheaper (depending on the approach taken) since it could require fewer trigonometric calculations than your approach. – Steven Stadnicki Sep 12 '13 at 6:33