This is an optimisation question rather than something I am unable to do. I have a naive, primitive working version of the code (I submit it below) but I am trying to make it work faster. I am asking for technical improvements.
For those familiar with Markov Chain Monte Carlo, this is an application of a random-walk Metropolis-Hastings step. For those not, the matter at hand is to apply the following step for all elements of set V for a function F(x) .
If F(proposed value) / F(existing value) < U, reject the proposed value (keep existing).
If F(proposed value) / F(existing value) > U, accept the proposed value (replace existing with proposed in the set).
or simply, in code form,
If[[F[Proposed]/[F[Existing] > U, V = Vnew, V = V]
where U is a drawing from a Uniform(0,1) distribution and the proposed value comes from a Standard Normal distribution. In the end, some elements will have been updated by proposed values while others will have remained the same.
The catch is that F(x) is a complicated expression that progresses over 5 sets V, Y, J, xV, xY of equal length T (so all variables have subscript t between 1 and T) . The proposed value is ONLY for $V_t$, and the function evaluates by taking "groups"
$F[(V_{t-1}, V_t, V_{t+1}), (Y_t, Y_{t+1}), (J_t, J_{t+1}), (xY_t, xY_{t+1}), (xV_t, xV_{t+1})]$
. In simple words, for evaluating at time t it considers both neighbouring values of V_t and the next values of $Y_t, J_t, xY_t, xV_t$.
Here is my version.
Omitting the ends for simplicity, I partitioned sets Y, J, xV, xY and then I formed a Do loop from 2 till T-1 where I partition V (Vfrac), I replace $V_t$ in set V with its proposed value and repartition (Vfrac2), I evaluate F(existing) / F (proposed) at time t (that is a simple [[i]] in the partitioned sets) and according to the result of the M-H step I replace (or not) SET V by the set that contains the proposed value. Please notice that the If operation returns SETS, not values.
Notice that the partitions take the form
(V1, V2, V3), V4, V5 ---> V1, (V2, V3, V4), V5 ---> V2, (V3, V4, V5),...
Y1, (Y2, Y3), Y4, Y5 ---> Y1, Y2, (Y3, Y4), Y5 ---> Y2, Y3, (Y4, Y5),... (all other sets)
as we cursor along the sets until T-1
This is dead slow and awkward, but it works. So I am open to any suggestion that will improve computational speed. Unfortunately I am not a good programmer, so I would appreciate some clarity, if that is not too much of a problem for you.
I supply a working version of the code where F is a toy function (if needed, I can supply the original). It might also be clearer than the description above.
Input
T := 100
Y = RandomVariate[NormalDistribution[0, 1], 100];
V := Table[1, {i, 1, T}];
J := RandomVariate[BernoulliDistribution[0.2], 100];
ξY := Table[-2.5, {i, 1, T}];
ξV := Table[1, {i, 1, T}];
μ := 0.05
ρ := -0.4
σV := 0.01
α := 0.02*0.73
β := -0.02
distV[{Y_, Y1_}, {ξY_, ξY1_}, {J_, JJ1_}, {ξV_, ξV1_}, {V0_, V_,
V1_}] :=
1/V Exp[-1/
2 ((V^2 - 2 (V0 + α + β*V0 + ξV*J) V -
2 ρ*Sqrt[σV] V (Y - μ - ξY*J))/(
σV (1 - ρ^2) V0) + (Y1 - μ - ξY1*JJ1)^2/
V + ((V1 - V - α - β*V - ξV1*JJ1) -
ρ*Sqrt[σV] (Y1 - μ - ξY1*JJ1))^2/(σV (1 - ρ^2) V))]
Loop construct - one Do loop inside anonther. The inner loop cursors along the length T sets from 2 till T, the outer loop repeats the process 250 times. V simply returns the output.
Do[
ξYξY = Partition[ξY, 2, 1, {-1, 1}, {}];
ξYfrac = Delete[ξYξY, 1];
JJ = Partition[J, 2, 1, {-1, 1}, {}];
Jfrac = Delete[JJ, 1];
ξVξV = Partition[ξV, 2, 1, {-1, 1}, {}];
ξVfrac = Delete[ξVξV, 1];
YY = Partition[Y, 2, 1, {-1, 1}, {}];
Yfrac = Delete[YY, 1];
{Do[{Venh = Append[Prepend[V, 0], 0];
Vfrac = Partition[Venh, 3, 1, {-3, 3}, {}];
Vprop = V[[i]] + RandomVariate[NormalDistribution[0, 0.05]];
If[Vprop > 0, Vprop = Vprop, Vprop = V[[i]]];
Vnew = ReplacePart[V, i -> Vprop];
Venh2 = Append[Prepend[Vnew, 0], 0];
Vfrac2 = Partition[Venh2, 3, 1, {-3, 3}, {}];
If[Log[
distV[Yfrac[[i]], ξYfrac[[i]], Jfrac[[i]], ξVfrac[[i]],
Vfrac2[[i]]]] -
Log[distV[Yfrac[[i]], ξYfrac[[i]], Jfrac[[i]], ξVfrac[[i]],
Vfrac[[i]]]] >
Log[RandomVariate[UniformDistribution[{0, 1}]]], V = Vnew,
V = V]}, {i, 2, T - 1}];
}, {250}];
V
Partitioning
you do at individual Y1, (Y2, Y3), Y4, Y5 ---> Y1, Y2, (Y3, Y4), Y5 make a function that takes the point and its both neighbours via VYJ[[All,i-1;;i+1]] ? $\endgroup$RandomVariate[UniformDistribution[{0, 1}]]
with justRandomReal[]
. $\endgroup$