# Generating particles from a distribution function using Monte Carlo

I have been given a 4D ($$x, y, v_x, v_y$$) distribution function, $$f(x,y,v_x, v_y)$$, generated by an external code. I want to generate a set of particles from this distribution function, say 10k particles each with their own $$x, y, v_x, v_y$$. I want to launch these 10k particles in my own code and follow them for some time in a particle pusher.

I hear some of you calling for a Monte Carlo method. People seem to specifically prefer the Metropolis-Hastings method for this kind of sampling. Here is my first question:

• In the Metropolis-Hastings method, if we reject a candidate we simply keep the previous sample. To my beginner eye this seems good if you want to reconstruct a distribution function, but not great if you want to generate unique samples from the distribution function. Is there a way to generate unique samples using the Metropolis-Hastings method? Or do I have to use some other method?

I have implemented a basic version of this method, based on an answer to a question here to first generate only $$x,y$$ --- I'll generate the 4D coordinates when I'm more confident:

import numpy as np
from scipy.stats import norm

# get spline for distribution function here
def MH(n = 10000, x_sigma = 0.5, y_sigma = 0.5):
x_cur, y_cur = 0.0, 0.0
x_innov = norm.rvs(size=n, scale=x_sigma)
y_innov = norm.rvs(size=n, scale=y_sigma)
u = np.random.uniform(size=n)
post_cur = spl(x_cur, y_cur)
posterior = np.zeros([n, 2])
accepted = 0
for t in range(n):
x_prop = x_cur + x_innov[t]
y_prop = y_cur + y_innov[t]
post_prop = spl(x_prop, y_prop)
alpha = post_prop / post_cur
if u[t] <= alpha:
x_cur = x_prop
y_cur = y_prop
post_cur = post_prop
accepted += 1
posterior[t,0] = x_cur
posterior[t,1] = y_cur
print('%.2f %% of points were accepted'%(accepted/n*100))
return posterior

samples = MH(n=10000)


Here spl is a spline of the input distribution function. This code seems to work well, generating something that looks like the original distribution function. A final question:

• I added another dimension by simply adding $$y$$, can I do the same with $$v_x, v_y$$ if I expand the spline into 4D?

Thanks! Any other related advice would be greatly appreciated.

• Is the density completely arbitrary, unbounded, and so forth? If there’s a bit of structure to it, you might explore an acceptance-rejection method. Mar 12 '21 at 16:54
• It's arbitrary, but always smooth and goes to zero at the boundaries too. Is there a particular acceptance-rejection method you think would be most suitable, or perhaps a favourite example? I take it Metropolis-Hastings isn't suitable?
– user38535
Mar 12 '21 at 16:56
• Out of curiosity is your $f$ built of one or several Gaussians? Simply guessing there’s a statistical mechanics connection. Mar 12 '21 at 17:56
• Please add a comment before downvoting a post. Mar 12 '21 at 19:27
• How complicated is the domain in which your parameters lie? Can it be meshed? Mar 12 '21 at 20:31

Simply redoing my earlier comment with a bit more space ...

Depending on the particulars, a simple acceptance-rejection method might be a good place to start. Suppose you know $$f$$ never gets bigger than $$f_{\text{max}}$$. Generate $$X$$ uniformly between $$x_{\text{min}}$$ and $$x_{\text{max}}$$, then similarly $$Y$$ from between $$y_{\text{min}}$$ and $$y_{\text{max}}$$ and so forth for $$V_x$$ and $$V_y$$. Compute $$F=f(X,Y,V_x,V_y)$$ and $$W$$ uniform between 0 and $$f_\text{max}$$. If $$W$$ is less than $$F$$ accept the point $$(X, Y, V_x, V_y)$$, otherwise repeat the process until you do. This generates points identically distributed as $$f$$ and independent, and might be what you need for your model.

If all you need is to compute the expectation of some quantities with respect to $$f$$, then you're doing Monte Carlo integration and other techniques such as importance sampling come into play, but I imagine you're on the right track already.

Here is a commonly used alternative:

Let's assume that your probability density $$f(x,y,v_x,v_y)$$ lives in a domain $$\Omega$$ that is bounded and that you can subdivide into "cells" $$\Omega_i$$. In one of the comments, you mention that $$\Omega$$ is simply a 4-dimensional box, and then the $$\Omega_i$$ could simply be subdivision into a regular mesh consisting of $$I=I_xI_yI_{v_x}I_{v_y}$$ small 4-dimensional sub-boxes.

Now compute the "relative probabilities" $$\tilde p_i=\int_{\Omega_i} f(x,y,v_x,v_y) \,dv_y\, dv_x\, dy\,dx$$. If $$f$$ is already a normalized probability distribution, then set $$p_i=\tilde p_i$$, otherwise $$p_i=\tilde p_i / \left(\sum_i \tilde p_i\right)$$. These $$p_i$$ correspond to the probability of finding a particle drawn from the original probability distribution in the $$i$$th cell. In practice, for most distributions $$f$$, we cannot compute $$\tilde p_i$$ exactly; but it can be approximated using the midpoint or trapezoidal rule and if the boxes $$\Omega_i$$ are small compared to the typical length scale of variation of $$f$$, then this approximation will be small.

Now, if you want to draw the $$j=1$$ (i.e., the first) particle from this distribution $$f$$, what you will do is draw a random number $$w_j$$ between zero and one, and you find that index $$k_j$$ so that $$\sum_{i=1}^{k_j-1} p_i \le w_j < \sum_{i=1}^{k_j} p_i$$, and that will be the cell in which the $$j$$th particle will be located. We will make the assumption that $$f$$ can be well approximated by a piecewise constant function, and draw the position $$\mathbf x_j$$ of the $$j$$th particle uniformly within the cell $$\Omega_{k_j}$$. Then repeat until you have as many particles as you wanted.

• Thanks for offering another method, it aids understanding of both methods suggested. If I understand correctly, the small 4D sub-box that has the relative probablity closed to the random number is selected? Presumably after ordering the boxes by their relative probability.
– user38535
Mar 15 '21 at 8:34
• Also, if this method was used to generate a sufficiently large number of particles, wouldn't some of these markers have identical 4D coordinates if the relative probabilities are sufficiently peaked and the grid resolution is small?
– user38535
Mar 15 '21 at 8:38
• Think about this in a 1d situation. You split the range in which your pdf is nonzero into small intervals. The $p_i$ correspond to the probability that a particle randomly from $f$ is in the $i$th interval. So the particle generation process first, for each particle, determines in which interval a particle should lie in, and then places it at a random location in this interval with a uniform probability. This is equivalent to approximating $f$ by a piecewise constant function. Mar 15 '21 at 9:54
• This is implemented here, by the way: github.com/dealii/dealii/blob/master/source/particles/… The documentation of that function is here: dealii.org/developer/doxygen/deal.II/… Mar 15 '21 at 10:09
• @PartPushed, Wolfgang's suggestion is a good one, probably the best. It's treating the density as a mixture model. In the literature it's referred to as a composition technique. Mar 15 '21 at 15:09