# Sampling a distribution (from a galaxy model)

I am reading the following article: http://www.kof.zcu.cz/st/dis/schwarzmeier/galaxy_models.html and am currently at section 5.6 (positions of bodies in a galaxy).

I am trying to redo the simulations myself in Python, but I have a questions regarding sampling of the distribution.

Given a distribution function (Hernquist distribution):

$$\rho(r)=\dfrac{M}{2\pi}\cdot \dfrac{a}{r(a+r)^3}$$

the article states that to simulate the distribution, one has to calculate the mass within a circle of radius $r$ like follows:

$$m(r) = \int_{0}^{r} 4 \pi r'^2 \rho(r') dr'$$

which is the cumulative mass distribution function.

The article states that this formula represents the PDF. However, looking at the shape this appears to me to be a CDF. Approaching infinity the function approaches 1.0 which is for me a clear indication of a CDF.

To sample this distribution the article cites the Von Neumann method where one has to generate an $r$ and a $m$ value, and scale them accordingly, and check whether or not they fall below the $m(r)$ graph. If they do they are accepted, else they are rejected.

Am I completely off by thinking this is wrong? If I do this I end up with the majority of stars ending up at the higher radii.

I have the feeling I am sampling an CDF here instead of a PDF. To get accurate results (e.g.: having the majority of the stars in the center) means that I have to perform the Von Neumann method with the $\rho(r)$ function.

I am unable to contact the author of the article, so that's why I am asking here.

The article looks indeed wrong. In fact, there are two mistakes.

First, you're right that the acceptance-rejection method has to be applied to $\rho(r)$, and not to $m(r)$. To understand how this idea works, suppose we want to generate a one-dimensional normalized distribution function $p(y)$. Now, let's assume we can rewrite this distribution function in terms of a variable $x$, such that it takes the form of a uniform distribution. That is, $$p(x) = \begin{cases} 1& \text{for 0\leqslant x \leqslant 1,}\\ 0& \text{elsewhere}. \end{cases}$$ Given $p(y)$, what is $x$? We have the Jacobian transformation $$p(y)dy = p(x(y))\left|\frac{dx}{dy}\right|dy = \left|\frac{dx}{dy}\right|dy,$$ which implies $$p(y) = \frac{dx}{dy},$$ assuming that $x(y)$ is an increasing function. Thus $$x = \int_0^y p(y')dy' = F(y).$$ In other words, the integral of $p(y)$ (or equivalently, the area under the curve) follows a uniform distribution. With this in mind, there are essentially two ways to perform a Monte-Carlo simulation.

The first way is the acceptance-rejection method: plot the curve $p(y)$ and uniformly generate a pair of numbers $(a,b)$ in the interval $([0,y_\max],[0,p_\max])$, where $y_\max$ and $p_\max$ are the upper bounds of $y$ and $p(y)$. If the coordinate $(a,b)$ lies under the curve $p(y)$, accept it; otherwise, reject it. If the coordinate is accepted, $y=a$ is generated point.

There are major drawbacks to this method: $y_\max$ and $p_\max$ can be infinite, so one would need a cut-off. And if $p(y)$ has a sharp peak, one ends up rejecting a lot of points.

A far more efficient method is to uniformly generate $x$, and calculate the corresponding $y$ by inverting $x=F(y)$: $$y = F^{-1}(x).$$ This automatically fills up the area under the curve, without rejecting points.

If the calculation of $F^{-1}(x)$ is too numerically involved, one can use a combination of both methods: introduce another (simpler) function $f(y)$ that lies everywhere above $p(y)$. Apply the inversion method to $f(y)$, generating a point $y$. Then uniformly generate a value $b$ in the interval $[0,f(y)]$. If $b\leqslant p(y)$, accept $y$; otherwise, reject it.

Now, consider the Hernquist distribution. Since it has a cusp at the origin, and the cumulative mass $m(r)$ is a simple function $$m(r) = M\frac{r^2}{(a+r)^2},$$ I'd definitely recommend the inversion method. But there is an important caveat here, and that's the second mistake in the article: $\rho(r)$ is not really a one-dimensional distribution. Instead it is a distribution in 3-dimensional space, and it is only a function of one variable due to spherical symmetry. In order to apply the Monte-Carlo method, we have to express $\rho$ as a truly one-dimensional distribution function, which we can do by expressing it in terms of the volume $$y = \frac{4\pi}{3}r^3.$$ Now we have $$p(y) = \rho(y) = \frac{M}{2\pi}\frac{a\,(3y/4\pi)^{-1/3}}{\left[a + (3y/4\pi)^{1/3}\right]^3},\\ F(y) = m(y) = \int_0^y\rho(y')dy' = M\frac{(3y/4\pi)^{2/3}}{\left[a + (3y/4\pi)^{1/3}\right]^2}.$$ Once we generated a point $y$, the corresponding radius is $$r=\left(\frac{3y}{4\pi}\right)^{1/3}.$$

There is an important consequence: there are likely more particles at large radii than around the centre, even though $\rho(r)$ is much larger at small radii. The reason is that particles between two radii $r$ and $(r+\Delta r)$ occupy a shell with volume $$V = \frac{4\pi}{3}\left[(r+\Delta r)^3-r^3\right].$$ The larger the radius $r$, the larger the volume of the shell, which means you need more particles to fill it and get $\rho(r)$. This is obvious in the case of a constant density, but it is also true for general densities.

• Thanks for the incredibly detailed answer! And I am glad I am right with the acceptance-rejection method :) – avanwieringen Jan 13 '14 at 7:07
• If I'm not mistaken, these images are from the Numerical Recipes text, which is copyrighted material, so you need to include references to where they are from. – Kyle Kanos Jan 28 '17 at 21:51

The CDF, $F(x)$, is related to the PDF, $f(x)$, via the relation:

$$F(x) = \int_{-\infty}^xdx'\,f(x')$$ In the case of radial distributions, your lower limit is obviously 0 and not $-\infty$. Thus, your CDF is $m(r)$ and the PDF is $4\pi(r')^2\rho(r')$ (technically should be $\rho(r')$ with the $4\pi(r')^2$ coming from $dx\to dr$ and spatial isotropy, but whatever).

In the case of invertible functions (e.g., $f(x)=Ax$ with normalization constant $A$), you can solve this as $$F(x)=\frac{A}{2}x^2\to x=\sqrt{\frac{2F}{A}}$$ If you generate a random number number, set this equal to $F$ and out pops your $x$ that satisfies the PDF.

In the case your functions is not invertible (often the case involving radial distributions), I'd suggest using a Newton-Raphson iterative method to finding roots for this case. This can be easily done since you know $F(x)$ (your CDF) $F'(x)=f(x)$ (your PDF): $$x_{new}=x_{old}-\frac{F(x_{old})}{f(x_{old})}$$ When $|x_{new}|<\epsilon$, then $x_{old}$ is the (approximate) root. Often, only a few iterations are necessary for convergence.

As an aside, I cannot stress enough how terrible a suggestion the acceptance-rejection method is. DO NOT USE THIS METHOD. It literally wastes computation time, something that ought to be viewed as precious, regardless of the prevalence of computers. Do not listen to anyone who says you have to use this method, they are completely and utterly wrong.

The Newton method I suggest above does not reject a single random number, will fit the PDF accurately, and is highly efficient.