2
$\begingroup$

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,

$$\frac{post(test)*prop(c_j)}{post(c_j)*prop(test)}>U(0,1)$$

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.

$\endgroup$

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.