Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

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}];