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$ Commented Mar 12, 2021 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$
    – user38535
    Commented Mar 12, 2021 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$ Commented Mar 12, 2021 at 17:56
  • 1
    $\begingroup$ Please add a comment before downvoting a post. $\endgroup$
    – boyfarrell
    Commented Mar 12, 2021 at 19:27
  • $\begingroup$ How complicated is the domain in which your parameters lie? Can it be meshed? $\endgroup$ Commented Mar 12, 2021 at 20:31

2 Answers 2


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$
    – user38535
    Commented Mar 15, 2021 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$
    – user38535
    Commented Mar 15, 2021 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$ Commented Mar 15, 2021 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$ Commented Mar 15, 2021 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$ Commented Mar 15, 2021 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.