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.