Difficulty in parallelizing function using package (MCMC)

Question

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

Answer 1

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}, ErrorBarPlots`ErrorBar @@ {#3, #4}} & @@@ 
  dat, PlotRangePadding -> {Scaled[0.15], Automatic}]

Indeed the Maximum likelihood method is unbiased

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