Fold/ NestList/ etc in MCMC loop for optimisation - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-09-22T04:56:06Z https://mathematica.stackexchange.com/feeds/question/186821 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/186821 2 Fold/ NestList/ etc in MCMC loop for optimisation Titus https://mathematica.stackexchange.com/users/30703 2018-11-27T23:33:08Z 2018-11-27T23:33:08Z <p>I am trying to simplify and optimise an old loop code I wrote. Since Fold and NestWhile use the results of the last iteration in the new iteration, I would like to use them in order to get rid of Do constructs in my version. The catch is certain operations. The example is a very simple template of a much more complicated structure.</p> <p>For those familiar with Bayesian inference/ MCMC codes, this mimics a hybrid Gibbs and Metropolis-Hastings sampler WITHOUT the statistical properties - I am only interested in the mechanics. For those who are not, I include a verbal description and a working coding example.</p> <p>I use SetDelayed to initialise parameters and vectors (variables) because they are used to store the value of the last iteration. I also define a posterior (target) and a proposal distribution, which will suggest candidate values. I define U to be used in a later control step and empty vectors for storing each loop's values. For simplicity the length of c and the number of iterations are set to 10. I know that some definitions are redundant, they help me track the process.</p> <pre><code>a := 1 b := 1 c := Table[1, {i, 1, 10}]; U := RandomVariate[UniformDistribution[{0, 1}]] astore := {} bstore := {} f[a_, b_, c_] := RandomVariate[NormalDistribution[a + Mean[c], Abs[b]]] g[a_, b_, c_] := RandomVariate[NormalDistribution[b + Mean[c], Abs[a]]] post[a_, b_, c_] := PDF[StudentTDistribution[Abs[b]], c] prop[a_, b_, c_] := PDF[NormalDistribution[a, Abs[b]], c] </code></pre> <p>The first two steps of each iteration of the loop are a drawing of a and b from distributions f and g. The third step is a cursoring over vector c (written as a loop-within-a-loop). For each element j in c a candidate (test) value is drawn and when the following control step is applied,</p> <p><span class="math-container">$$\frac{post(test)*prop(c_j)}{post(c_j)*prop(test)}&gt;U(0,1)$$</span></p> <p>If the test holds, then <span class="math-container">$c_j$</span> is replaced by test in vector c. If not, nothing happens (logarithms are for numerical reasons).</p> <pre><code>Do[{a = f[a, b, c], b = g[a, b, c], Do[{ex = c[[j]], test = RandomVariate[NormalDistribution[a, Abs[b]]], If[Log[post[a, b, test]*prop[a, b, ex]] - Log[post[a, b, ex]*prop[a, b, test]] &gt; Log[U], c = ReplacePart[c, j -&gt; test], Nothing]}, {j, 1, 10}], astore = Append[astore, a], bstore = Append[bstore, b] }, {i, 1, 10}] </code></pre> <p>The code must store a and b for each iteration and update vector c properly. </p> <p>I would like to 1) Improve on speed (so this can be a horse-race between different constructs) 2) Use Fold and/ or NestList which are made for such operations.</p> <p>Any other improvement is also welcome. Finally, there is a considerable complication in c that includes walkers, but for now I will stick to that.</p> <p>If you require any further information by all means, ask and I will add it. </p>