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