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.