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.

  • $\begingroup$ 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. $\endgroup$ – A rural reader Mar 12 at 16:54
  • $\begingroup$ 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? $\endgroup$ – PartPushed Mar 12 at 16:56
  • $\begingroup$ Out of curiosity is your $f$ built of one or several Gaussians? Simply guessing there’s a statistical mechanics connection. $\endgroup$ – A rural reader Mar 12 at 17:56
  • 1
    $\begingroup$ Please add a comment before downvoting a post. $\endgroup$ – boyfarrell Mar 12 at 19:27
  • $\begingroup$ How complicated is the domain in which your parameters lie? Can it be meshed? $\endgroup$ – Wolfgang Bangerth Mar 12 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.

  • $\begingroup$ 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. $\endgroup$ – PartPushed Mar 15 at 8:34
  • $\begingroup$ 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? $\endgroup$ – PartPushed Mar 15 at 8:38
  • $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Mar 15 at 9:54
  • $\begingroup$ 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/… $\endgroup$ – Wolfgang Bangerth Mar 15 at 10:09
  • $\begingroup$ @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. $\endgroup$ – A rural reader Mar 15 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.