Monte Carlo code optimization

• To: mathgroup at smc.vnet.net
• Subject: [mg130556] Monte Carlo code optimization
• From: mccarter.patrick at gmail.com
• Date: Sat, 20 Apr 2013 05:47:10 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• Delivered-to: l-mathgroup@wolfram.com
• Delivered-to: mathgroup-newout@smc.vnet.net
• Delivered-to: mathgroup-newsend@smc.vnet.net

Hi All,

I am working on an optimization problem using a Markov Chain Monte Carlo approach. I have a crude for mathematica implementation using for loops. However I would really like to learn how to:
1) transition from a procedural to functional approach to implementing the algorithm. Can anyone offer hints at how to optimize the notebook towards a more functional approach.
2) I'd also like to add a time dependent signal to the NDSolve function. However I am not quite sure how to do this, if I use the If statement, then I must add the If to each equation. This seems inefficient. Is there a bette r approach to adding time dependent factors to the ODE solver. Say I want to wait until steady state to add the signal?

Thanks,
I have attached the notebook with comments. If any more are needed please let me know. Thanks for your assistance and pointers!

Patrick

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;

(* I'd like to turn on a stimulus for a set amount of time *)
(* \
If[25<t<50, stim=0.01, Else stim = 1.0; *)

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[bestParamSets]
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, range}];
kb[[0]] = RandomReal[{0.0, range}];
decayscore[[0]] = chi2[ka[[0]], kb[[0]]]
beta = -(Log[0.2]/decayscore[[0]])*10^2

For[i = 1., i < N[runlength], i++,
{(* Open Main For Loop *)

ka[[i]] =
ka[[i - 1]] + weight*ka[[i - 1]]*RandomReal[{-1.0, 1.0}]; (*
Other distributions are possible *)

kb[[i]] =
kb[[i - 1]] + weight*kb[[i - 1]]*RandomReal[{-1.0, 1.0}];
(*************************************************)
(*
This is where the evaluation of the parameter sets begins *)
(*************************************************)

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[-(1/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 < v *)
{
ka[[i + 1]] = ka[[i]],
kb[[i + 1]] = kb[[i]],
AppendTo[
bestParamSets, {decayscore[[i]], v[[i]], ka[[i]], kb[[i]],
i}];
},
{
ka[[i]] = ka[[i - 1]],
kb[[i]] = kb[[i - 1]],
AppendTo[
thrownsets, {decayscore[[i]], v[[i]], ka[[i]], kb[[i]], i}];
}
]
} (* Close Main For Loop // How long does this loop take to run? *)
\
] // AbsoluteTiming

Length[bestParamSets]
Length[thrownsets]

N[Length[bestParamSets]/(Length[bestParamSets] + Length[thrownsets])]

bps = Ordering[bestParamSets, sets];
bestsets = Table[bestParamSets[[j]], {j, bps}];
kA = Flatten[Table[bestParamSets[[i, 3]], {i, 1, Length[bestsets]}]];
kB = Flatten[Table[bestParamSets[[i, {4}]], {i, 1, Length[bestsets]}]];
pSets = Table[bestsets[[i, {3, 4}]], {i, 1, sets}];

• Prev by Date: Re: Easy to get the audio out of sync with the graphics (Repost)
• Next by Date: Re: how to define and analyze function with multiple
• Previous by thread: Re: Find the function that best fits
• Next by thread: Re: Parametric List Plot??