Generating particles from a distribution function using Monte Carlo

Question

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.

Answer 1

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.

Answer 2

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.