# Fast implementation of Metropolis-Hastings update with multiple sets

This is an optimisation question rather than something I am unable to do. I have a naive, primitive working version of the code (I submit it below) but I am trying to make it work faster. I am asking for technical improvements.

For those familiar with Markov Chain Monte Carlo, this is an application of a random-walk Metropolis-Hastings step. For those not, the matter at hand is to apply the following step for all elements of set V for a function F(x) .

If F(proposed value) / F(existing value) < U, reject the proposed value (keep existing).

If F(proposed value) / F(existing value) > U, accept the proposed value (replace existing with proposed in the set).

or simply, in code form,

If[[F[Proposed]/[F[Existing] > U, V = Vnew, V = V]

where U is a drawing from a Uniform(0,1) distribution and the proposed value comes from a Standard Normal distribution. In the end, some elements will have been updated by proposed values while others will have remained the same.

The catch is that F(x) is a complicated expression that progresses over 5 sets V, Y, J, xV, xY of equal length T (so all variables have subscript t between 1 and T) . The proposed value is ONLY for $V_t$, and the function evaluates by taking "groups"

$F[(V_{t-1}, V_t, V_{t+1}), (Y_t, Y_{t+1}), (J_t, J_{t+1}), (xY_t, xY_{t+1}), (xV_t, xV_{t+1})]$

. In simple words, for evaluating at time t it considers both neighbouring values of V_t and the next values of $Y_t, J_t, xY_t, xV_t$.

Here is my version.

Omitting the ends for simplicity, I partitioned sets Y, J, xV, xY and then I formed a Do loop from 2 till T-1 where I partition V (Vfrac), I replace $V_t$ in set V with its proposed value and repartition (Vfrac2), I evaluate F(existing) / F (proposed) at time t (that is a simple [[i]] in the partitioned sets) and according to the result of the M-H step I replace (or not) SET V by the set that contains the proposed value. Please notice that the If operation returns SETS, not values.

Notice that the partitions take the form

(V1, V2, V3), V4, V5 ---> V1, (V2, V3, V4), V5 ---> V2, (V3, V4, V5),...

Y1, (Y2, Y3), Y4, Y5 ---> Y1, Y2, (Y3, Y4), Y5 ---> Y2, Y3, (Y4, Y5),... (all other sets)

as we cursor along the sets until T-1

This is dead slow and awkward, but it works. So I am open to any suggestion that will improve computational speed. Unfortunately I am not a good programmer, so I would appreciate some clarity, if that is not too much of a problem for you.

I supply a working version of the code where F is a toy function (if needed, I can supply the original). It might also be clearer than the description above.

Input

T := 100
Y = RandomVariate[NormalDistribution[0, 1], 100];
V := Table[1, {i, 1, T}];
J := RandomVariate[BernoulliDistribution[0.2], 100];
ξY := Table[-2.5, {i, 1, T}];
ξV := Table[1, {i, 1, T}];

μ := 0.05
ρ := -0.4
σV := 0.01
α := 0.02*0.73
β := -0.02

distV[{Y_, Y1_}, {ξY_, ξY1_}, {J_, JJ1_}, {ξV_, ξV1_}, {V0_, V_,
V1_}] :=
1/V Exp[-1/
2 ((V^2 - 2 (V0 + α + β*V0 + ξV*J) V -
2 ρ*Sqrt[σV] V (Y - μ - ξY*J))/(
σV (1 - ρ^2) V0) + (Y1 - μ - ξY1*JJ1)^2/
V + ((V1 - V - α - β*V - ξV1*JJ1) -
ρ*Sqrt[σV] (Y1 - μ - ξY1*JJ1))^2/(σV (1 - ρ^2) V))]


Loop construct - one Do loop inside anonther. The inner loop cursors along the length T sets from 2 till T, the outer loop repeats the process 250 times. V simply returns the output.

Do[

ξYξY = Partition[ξY, 2, 1, {-1, 1}, {}];
ξYfrac = Delete[ξYξY, 1];
JJ = Partition[J, 2, 1, {-1, 1}, {}];
Jfrac = Delete[JJ, 1];
ξVξV = Partition[ξV, 2, 1, {-1, 1}, {}];
ξVfrac = Delete[ξVξV, 1];
YY = Partition[Y, 2, 1, {-1, 1}, {}];
Yfrac = Delete[YY, 1];

{Do[{Venh = Append[Prepend[V, 0], 0];
Vfrac = Partition[Venh, 3, 1, {-3, 3}, {}];
Vprop = V[[i]] + RandomVariate[NormalDistribution[0, 0.05]];
If[Vprop > 0, Vprop = Vprop, Vprop = V[[i]]];
Vnew = ReplacePart[V, i -> Vprop];
Venh2 = Append[Prepend[Vnew, 0], 0];
Vfrac2 = Partition[Venh2, 3, 1, {-3, 3}, {}];
If[Log[
distV[Yfrac[[i]], ξYfrac[[i]], Jfrac[[i]], ξVfrac[[i]],
Vfrac2[[i]]]] -
Log[distV[Yfrac[[i]], ξYfrac[[i]], Jfrac[[i]], ξVfrac[[i]],
Vfrac[[i]]]] >
Log[RandomVariate[UniformDistribution[{0, 1}]]], V = Vnew,
V = V]}, {i, 2, T - 1}];

}, {250}];
V

• Titus, two suggestions. (1) you are applying a lot of magic with Partition, Append and Prepend and in the end it results into a complex datastructure that is being manipulated in a way that is difficult (at least for me, but also for a Compiler) to understand. If you can simplify this, it will make a lot of manipulations obsolete and open the door to speeding this thing up, (2) please look into the differences between := and =, the former is SetDelayed, which is a.o. used to define Functions. It is no use defining constants this way. Apr 20, 2016 at 10:56
• Related: 111224.
– gwr
Apr 20, 2016 at 15:39
• Would it be possible to construct a 3x100 array of {{V1,Y1,J1},{{V2,Y2,J2},... ,{{V100,Y100,J100}} - let's call it VYJ - and instead of the complex Partitioning you do at individual Y1, (Y2, Y3), Y4, Y5 ---> Y1, Y2, (Y3, Y4), Y5 make a function that takes the point and its both neighbours via VYJ[[All,i-1;;i+1]] ? Apr 22, 2016 at 1:25
• Not a great improve in speed, but you can replace RandomVariate[UniformDistribution[{0, 1}]] with just RandomReal[]. Apr 28, 2017 at 1:09

Append and Prepend and possibly ReplacePartare most likely slowing your code down substantially. I could not recode your stuff to work without these constructs withing the limited time I had available, nor do I understand competely why you approach the issue this way. Regardless, there is a good and clear demonstration by Philip Gregory that gives an excellent example implementation of MCMC.

"Markov Chain Monte Carlo Simulation Using the Metropolis Algorithm" from the Wolfram Demonstrations Project http://demonstrations.wolfram.com/MarkovChainMonteCarloSimulationUsingTheMetropolisAlgorithm/

Contributed by: Philip Gregory (Physics and Astronomy, University of British Columbia)

I am sure this is a clear enough example for you.

Try to recode your example in this way, e.g. replace Append, Prepend and ReplacePart and then Compile your code to get another performance boost.

To demonstrate performance improvement, I am stealing this example, and remove the manipulate (there's no need to explain this construct for your question) and I hope that if what I do below, you can apply to your code as follows, you will see a substantial improvement as well:

1. Work with Nest or NestList and get rid of Append, Prepend and ReplacePart;
2. Compile;
3. Compile to C.

We then get something like this for the uncompiled version:

mcmcfun = Function[s,
Module[{xm, ym, proposal, xp, yp, p2, p1, proposalSigma = 0.2,
pdf = 0.087 E^(
1/2 (-(x - 4) (0.595 (x - 4) -
0.238 (y - 3)) - (-0.238 (x - 4) + 0.595 (y - 3))  (y -
3))) + E^(1/2 (-x^2 - y^2))/(2 \[Pi])},
{xm, ym} = s;
xp = RandomReal[NormalDistribution[xm, proposalSigma]];
yp = RandomReal[NormalDistribution[ym, proposalSigma]];
p2 = pdf /. {x -> xp, y -> yp};
p1 = pdf /. {x -> xm, y -> ym};
proposal = p2/p1;
If[RandomReal[] <= proposal
, {xp, yp}
, {xm, ym}]]];


And this is called 1 Million times by:

(sim = NestList[
mcmcfun[#] &,
{-4,9},
1000000];) // AbsoluteTiming

{132.823563, Null}


It takes roughly 2 minutes on an Intel i7 Macbook pro. Depending on your purpose, this may be acceptable performance. If not, continue reading.

Speeding it up 60 times by compiling it and 240 times by compiling to C by using the option CompilationTarget -> "C" as follows:

pdf = Compile[{{x, _Real}, {y, _Real}},
0.087 E^(
1/2 (-(x - 4) (0.595 (x - 4) -
0.238 (y - 3)) - (-0.238 (x - 4) + 0.595 (y - 3))  (y -
3))) + E^(1/2 (-x^2 - y^2))/(2 \[Pi]),
Parallelization -> True, CompilationTarget -> "C"];

mcmcfunComp = Compile[{{s, _Real, 1}},
Module[{xm, ym, proposal, xp, yp, p2, p1, proposalSigma = 0.2},
{xm, ym} = s;
xp = RandomReal[NormalDistribution[xm, proposalSigma]];
yp = RandomReal[NormalDistribution[ym, proposalSigma]];
p2 = pdf[xp, yp];
p1 = pdf[xm, ym];
proposal = p2/p1;
If[RandomReal[] <= proposal
, {xp, yp}
, {xm, ym}]]
, CompilationOptions -> {"InlineExternalDefinitions" -> True},
Parallelization -> True, CompilationTarget -> "C"];


This takes less than 0.5 seconds!

(sim = NestList[
mcmcfunComp[#] &,
{-4,9},
1000000];) // AbsoluteTiming

{0.464945, Null}


Please note, I included a Parallelization -> True option, this does not speed things up as there is nothing to parallelise. However I assume you will want to use several walkers; in that case the parallelization will likely help you speed up things further.

Hope this helps.

• This gave me a nice idea to get rid of the uncomfortable Do loop and write the whole thing in a more compact form. However, I am a bit confused about the Nest/ NestList function you suggest. NestList applies function F repeatedly over each element of list x, while I want it to be applied only once. I could probably use Map or MapAt (my version relies on their intuition actually) but I am more interested in your solution. Apr 18, 2016 at 8:00
• @Titus That's good to hear. I am not sure but it seems you refer to FoldList when you say NestList. I am using NestList to feed the Function mcmcfun the initial state {-4,9} , and the result of that, I feed into mcmcfun again, and the result of that again ..... and I do this 1 million times.... If however I wanted to feed a list, I would indeed use FoldList. Let me know if I understood you correctly. Apr 19, 2016 at 0:36
• @Titus Perhaps I have different understanding of MCMC; The use of Map would work well for Monte Carlo simulations (although there are faster ways), but the "goal seeking" nature of MCMC is IMHO better served by using NestList, as it feeds the results back into the mcmcfun, where Map would require an additional construct to do this. Apr 19, 2016 at 1:08
• If I am reading this correctly, NestList would be suitable if I wanted to draw a single value e.g. from Gibbs sampling. That way NestList would indeed feed the results back into mcmcfun and also deliver all realisations. However, in my case, the entire cursoring over set V takes place in one step (walking over triads V_t-1 , V_t , V_t+1 as I explain in my initial post leads (see the partition example) to a different set V that is ONE drawing for my purposes). I will check to see how MapAt would work with your example. Does that clarify anything? Apr 19, 2016 at 10:14
• NestList is not restricted to a single value, you can feed - and actually I think this is what you should do - the function your 5 sets, and rewrite your function so you will get the new sampled sets as a result so you can feed it again as the next iteration. It will help me if you can provide a minimal working example, including input and output so I can update my answer. Apr 19, 2016 at 12:59