For the following code, which runs independent chains of Random Walk Metropolis, there is the obvious potential for parallelization. I am wondering what strategies people would recommend. This is related to a prior question I had about how to use pmap
when it was being applied to something which was an in place transformation.
function Boltzmann_likelihood(x, V, beta)
w = exp(-beta * V(x));
return w
end
function RWM!(X0, V::Vf, beta; niterations=10^2, Dt = 1.e-1) where {Vf}
Xp = similar(X0);
naccept = 0;
p0 = Boltzmann_likelihood(X0, V, beta);
gaussian_coef = sqrt(2 * Dt/beta);
for j = 1:niterations
@. Xp = X0 + gaussian_coef * randn();
pp = Boltzmann_likelihood(Xp, V, beta);
a = min(1, pp/p0);
#a = 1.
if rand()<a
naccept = naccept+1;
@. X0 = Xp;
p0 = pp;
end
end
X0
end
nsamples = 10^2;
d= 2;
beta = 1;
function V(x)
return (dot(x,x)-1)^2
end
X0vals = [randn(d) for i=1:nsamples];
Xbar= zeros(d);
for k = 1:nsamples
RWM!(X0vals[k], V, beta)
@. Xbar += X0vals[k]/nsamples;
end
println(Xbar);