Hello All,
I am on a journey to become adept at developing functional programs with Mathematica. I have a short prototype
code that I am using to learn how to convert procedural programs into functional programs and am wondering if anyone
could provide any more advise? Specifically I would like to learn how to write this code in a much more memory friendly manner.
And also learn a better way to run the main loop of the Markov Chain Monte Carlo simulation that minimizes the distance between
the experimental radioactive decay data and my model. Does anyone have any suggestions on this or how to develop into better functional programmer.
Thanks,
Pat
Clear["Global`*"]
texp = {0.2, 2.2, 4.0, 5.0, 6.0, 8.0, 11.0, 15.0, 18.0, 26.0, 33.0,
39.0, 45.0};
dexp = {35., 25., 22.1, 17.9, 16.8, 13.7, 12.4, 7.5, 4.9, 4.0, 2.4,
1.4, 1.1};
dt = Transpose[{texp, dexp}];
dataPlot = ListPlot[
dt,
PlotStyle -> {Red, AbsolutePointSize[6]},
Frame -> {True, True, False, False}, Axes -> False,
FrameLabel -> {"Time", "Species A"},
LabelStyle -> {FontSize -> 20, FormatType -> Bold},
RotateLabel -> False,
ImageSize -> {550}];
runlength = 10000;
sets = runlength/10.;
testing = 10.0;
Mean[dexp];
StandardDeviation[dexp];
stim = 0.01;
chi2[ka_?NumberQ, kb_?NumberQ] := Block[{sol, A, B, c}, sol =
NDSolve[
{ A'[t] == -ka*A[t],
B'[t] == ka*A[t] - kb*B[t],
c'[t] == kb*B[t],
A[0] == 35.,
B[0] == 0.,
c[0] == 0.},
{A, B, c},
{t, 0., 100.}][[1]];
Apply[Plus, (dexp - (A[t] /. sol /. t -> texp))^2]];
DateString[];
RandomSeed[AbsoluteTime[%]]; (* Set Random Seed based on current time *)
Clear[deltaE]
Clear[freeE]
Clear[ka]
Clear[kb]
Clear[u, v]
bestParamSets = {};
thrownsets = {};
(* These arrays are for storage *)
ka = ConstantArray[0, {runlength}];
kb = ConstantArray[0, {runlength}];
deltaE = ConstantArray[0, {runlength}];
freeE = ConstantArray[0, {runlength}];
v = ConstantArray[0, {runlength}];
u = ConstantArray[0, {runlength}];
decayscore = ConstantArray[0, {runlength}];
(**************************************************)
(******** Starting the Metropolis Algorithm ********)
(**************************************************)
range = 10.0;(* Pick different initial guesses for param. sets *)
weight = 1./10;
ka[[0]] = RandomReal[{0.0, 0.2}];
kb[[0]] = RandomReal[{0.0, 0.2}];
least = 0.0004;
max = 1.000;
decayscore[[0]] = chi2[ka[[0]], kb[[0]]];
beta = 1/.01;
(tbl = Reap[
Table[(*Open Main For Loop*)
ka[[i]] = ka[[i - 1]] + weight*ka[[i - 1]]*RandomReal[{-1.0, 1.0}];
kb[[i]] = kb[[i - 1]] + weight*kb[[i - 1]]*RandomReal[{-1.0, 1.0}];
(*************************************************)
(*This is where the evaluation of the parameter sets begins*)
(*************************************************)
If[
least < ka[[i]] && ka[[i]] < max &&
least < kb[[i]] && kb[[i]] < max,
decayscore[[i]] = chi2[ka[[i]], kb[[i]]];(*This is the function evaluated for each parameter set*)
deltaE[[i]] = decayscore[[i]] - decayscore[[i - 1]];(*The change in energy for each step*)
freeE[[i]] = Exp[-beta*deltaE[[i]]];(*This is a free energy associated with each step*)
v[[i]] = Min[1, freeE[[i]]];
u[[i]] = RandomReal[{0, 1}];
(**************************************************)
(******** Metropolis Selection Algorithm ********)
(**************************************************)
If[u[[i]] < v[[i]],(*If u<v always accept that parameter set, else accept with probability u {data}]*);
ka[[i + 1]] = ka[];
kb[[i + 1]] = kb[];
Sow[{True, decayscore[], v[], ka[], kb[], i}],(*bestParamSets*)
ka[] = ka[[i - 1]];
kb[] = kb[[i - 1]];
],
ka[] = ka[[i - 1]];
kb[] = kb[[i - 1]]
],
{i, 1,runlength - 1}]][[2, 1]];
bestParamSets = Cases[tbl, {True, data__} -> {data}];
) // AbsoluteTiming (*Close Main For Loop How long does this loop take to run?*)
