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.

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.

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.

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]

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,


If the test holds, then $c_j$ is replaced by test in vector c. If not, nothing happens (logarithms are for numerical reasons).

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]] > Log[U], 
 c = ReplacePart[c, j -> test], Nothing]}, {j, 1, 10}],

astore = Append[astore, a],
bstore = Append[bstore, b]

}, {i, 1, 10}]

The code must store a and b for each iteration and update vector c properly.

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.

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.

If you require any further information by all means, ask and I will add it.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.