2

I am trying to use compile to speed up my Gibbs sampler. I know I can use Reap and Sow to speed up tremendously. I have a more complicated Bayesian posterior that I want to speed up with compile so I thought id start with the simpler Bivariate sampler to use compile.

My original code is

targetMean = {0, 0};(*Mean of Target Distribution*)
covarianceParam = {.5, .5} ; (*Off diagnal terms of covariance matrix*)
dimensions = {1, 2}; (*Index of each dimension*)

paramSet1 = {
          nSamples -> 50000,
          m0 -> targetMean,
          r0 -> covarianceParam
         };

mCond[x_, d_] := (m0 /. paramSet1)[[d]] + (r0 /. paramSet1)[[d]]
                 *(x- (m0 /. paramSet1)[[Cases[{1, 2}, Except[d]][[1]]]]);
rCond[d_] := (Sqrt[1 - (r0 /. paramSet1)[[d]]^2]);
xCond[m_, r_] := RandomVariate[NormalDistribution[m, r]]


Gibbs[initx0_, ParamList0_] :=
    Module[
       {initx = initx0,
        intr = initr0,
        intTh = intTh0,
        paramSet = ParamList0,
        im, ir, ix}, (*Defines varibles for for function*)
        ix = initx; (*Another name for varibles shorter, can see what they mean*)

        xArray = {}; (*creates and empty array to store future data*)

    Monitor[For[t = 1, t < nSamples /. paramSet, t = t + 1 /. paramSet,

     For[iD = 1, iD < 3, iD++,

        im] = mCond[ix, iD] /. paramSet;
        ir = rCond[iD] /. paramSet;
        ix = xCond[im, ir] /. paramSet;


        AppendTo[xArray, ix]

       ];
     ], {t, ProgressIndicator[t, {1, nSamples}] /. paramSet}]
    ]

I then tried to use compile in the following code. I'm still not sure exactly how to compile (not the best documentation on the internet). I would love to get it to work for my Bayesian Gibbs Sampler since that takes quite a bit of time.

 m0 = {0, 0};(*Mean of Target Distribution*)
 r0 = {.5, .5} ; (*Off diagonal terms of covariance matrix*)
 dimensions = {1, 2}; (*Index of each dimension*)
 nSamples = 10000;

mCond[x_, d_] := m0[[d]] + r0[[d]]*(x - m0[[Cases[{1, 2}, Except[d]][[1]]]]);
rCond[d_] := (Sqrt[1 - r0[[d]]^2]);
xCond[m_, r_] := RandomVariate[NormalDistribution[m, r]];


Gibbs = Compile[{{ix1, _Real}},

Module[
    {ix = ix1,
     ir, im},

    xArray = {};

    Monitor[For[t = 1, t < nSamples, t = t + 1,

     Do[

        im = Evaluate[mCond[ix, iD]];
        ir = Evaluate[rCond[iD]];
        ix = Evaluate[xCond[im, ir]];

        AppendTo[xArray, ix];
        , {iD, 1, 2}];
     ], {t, ProgressIndicator[t, {1, nSamples}]}]
 ]]

I get the error

CompiledFunction::cfse: Compiled expression Null 
   should be a machine-size real number. >>

and

CompiledFunction::cfex: Could not complete external evaluation at instruction 3; 
  proceeding with uncompiled evaluation. >>
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
ADDonut
  • 139
  • 6
  • Here is a list of compilable functions in Mathematica. As you can see, none of Monitor, NormalDistribution or importantly Pattern are there. You can issue << CompiledFunctionTools\ in a notebook and useCompilePrinton your compiled attempt to see where it breaks of compile to do "normal" _Mathematica_ stuff; you will seeMainEvaluate` in the output. – Marius Ladegård Meyer Jun 14 '16 at 21:55
  • 2
    Note however that there are probably lots you can do with your code to speed it up before Compile becomes a relevant optimization. This is a good place to start; don't use AppendTo, don't use For loops, try to have your functions act on whole lists of data instead of element-wise operations etc. – Marius Ladegård Meyer Jun 14 '16 at 21:59
  • 2
    @Marius, actually the generation of normal variates is a compilable operation; see for instance CompilePrint[Compile[{{_Real}}, RandomVariate[NormalDistribution[], 2]]]. – J. M.'s missing motivation Jun 14 '16 at 23:34
  • 1
    @MariusLadegårdMeyer For more detail about the compilation, you may want to read this: http://mathematica.stackexchange.com/q/1124/1871 – xzczd Jun 15 '16 at 05:57
  • I got compile to work thats for all your help. Great resources. The next thing I need to figure out is how to get ride of my for loops. I think that could really speed things up. Any tips? – ADDonut Jun 15 '16 at 22:20

0 Answers0