# Difficulty in parallelizing function using package (MCMC)

## Difficulty in parallelizing Packages using DistributeDefinitions

Context

I am trying to parallelize some parameter estimation using MCMC.

As a test run I expect it to find the width of a Gaussian.

If I use the built in Mathematica function FindDistributionParameters, it parallelizes nicely:

Table[x = RandomVariate[NormalDistribution[0, 2/10], 50000];
FindDistributionParameters[x, NormalDistribution[0, s]], {i,
40}] // Timing


{0.324758, {{s -> 0.200035},…}

runs much slower than

ParallelTable[x = RandomVariate[NormalDistribution[0, 2/10], 50000];
FindDistributionParameters[x, NormalDistribution[0, s]], {i,
40}] // Timing


{0.016444, {{s -> 0.199899}, …}

If I now use a MCMC package (following this answer)

<<MCMC.m

Unprotect[Monitor];
Monitor = # &;
Protect[Monitor];
DistributeDefinitions[MCMC];


For some reason this (non parallel) implementation is as slow

Table[x = RandomVariate[NormalDistribution[0, 2/10], 500];
f[s_] = LogLikelihood[NormalDistribution[0, s], x];
MCMC[f[s], {{s, 0.1, 0.4, Reals}}, 1000]["BestFitParameters"]
, {i, 4}] // Timing


{2.26298,{{s->0.188522},{s->0.208759},{s->0.208809},{s->0.1983}}

as that (parallel) implementation:

ParallelTable[x = RandomVariate[NormalDistribution[0, 2/10], 500];
f[s_] = LogLikelihood[NormalDistribution[0, s], x];
MCMC[f[s], {{s, 0.1, 0.4, Reals}}, 1000]["BestFitParameters"]
, {i, 4}] // Timing


{3.38599,{{s->0.199691},{s->0.200183},{s->0.206191},{s->0.201176}}

## Question

Why does ParallelTable not scale up on this problem given that the different jobs are fully independent?

PS: note that an alternative option could involve parallelizing the function itself, following e.g. this answer

• I recently found that in using packages I would need to use Needs[] as well as ParallelNeeds[] to have the definitions distributed correctly to the Kernels. – gwr Mar 3 '17 at 15:58
• Suprised the heck out of me, too... just when one thought, one knew all the tricks. :) – gwr Mar 3 '17 at 16:02
• When I do this with MCMC, I get Min::nord and FrontEndObject::notavail (albeit getting the correct result). Do you get the same messages? – gwr Mar 3 '17 at 16:17
• @gwr yes unless I do ParallelDo[Unprotect[Monitor]; Monitor = # &; Protect[Monitor]; , {4}]; because the function uses Monitor. – chris Mar 3 '17 at 16:36
• I'm not an expert but "parallel" and "random numbers" ring some bells. I think there are some peculiarities with that combination. If unsure you might want to make sure that what you do is correct... – Albert Retey Mar 4 '17 at 20:19

## Loading Packages with ParallelNeeds

As Szabolcs has nicely posted here distributing definitions to the Kernels by using DistributeDefinitions like so:

DistributeDefinitions[ "PackageContext"]


will not work. Since definitions of the default context are distributed automatically, one needs to make sure, that packages loaded via Needs also "reach" the Kernels.

This can be achieved by using ParallelNeeds. But one still needs to load the package into the global context thus in your case:

Needs[ "MCMC"];
ParallelNeeds[ "MCMC"];


That should do the trick.

In practice for the problem at hand:

ParallelDo[Unprotect[Monitor];Monitor = # &; Protect[Monitor];, {4}]


Then we can compute best estimate and error on estimate as follow:

res = ParallelTable[
x = RandomVariate[NormalDistribution[0, 2/10], 5000];
f[m_, s_] = LogLikelihood[NormalDistribution[m, s], x];
MCMC[f[m, s], {{m, 0.1, 0.01, Reals},
{s, 0.1, 0.01, Reals}}, 5000][{"BestFitParameters",
"ParameterErrors"}]
, {i, 32*8}];

dat = Join[{m, s} /. Map[#[[1, 2]] &, res] //
Transpose, {m, s} /. Map[#[[2, 2]] &, res] // Transpose] //
Transpose;
Needs["ErrorBarPlots"];
ErrorListPlot[{{#1, #2}, ErrorBarPlotsErrorBar @@ {#3, #4}} & @@@
dat, PlotRangePadding -> {Scaled[0.15], Automatic}]


Indeed the Maximum likelihood method is unbiased

 dat1 = {m, s} /. Map[#[[1, 2]] &, res]; dat1 // DensityHistogram
`