For some time now, I have used this procedural programming approach to fit data in my research that is too cumbersome for the built in Mathematica functions. (Although this simple example works fine for the built in function, the complex problems that I am actually solving do not. More explicitly, those problems often fail to converge to a good solution, even with good starting conditions for the model parameters.) I now have some time to really think about my codes and would really appreciate some direction to make them much more mathematica efficient. I use a Markov Chain Monte Carlo algorithm to search for model parameters that minimize the difference between my model generated curves and the experimental data via a Chi squared scoring function. I've attached a short example of the types of problems that I solve with this approach written in the aforementioned procedural style. I've seen some instances of vectorizing mathematica code to make it much more efficient, however I am unfamiliar with how to do this myself. Thus far, I've tried this by calling all of the Random numbers before the crux of the algorithm, however, I did not achieve any speedup. At this point, I am a bit miffed as to what went wrong, or what I could try next. The main algorithm starts in the section prefaced by the Starting Metropolis Algorithm in the code. I've added the entire notebook for the sake of completeness. if anyone has any suggestions or expertise I would be greatly appreciative!


texp = {0.2, 2.2, 4., 5., 6., 8., 11., 15., 18., 26., 33., 39., 45.};  
dexp = {35., 25., 22.1, 17.9, 16.8, 13.7, 12.4, 7.5, 4.9, 4., 2.4, 1.4, 1.1};  
dt = Transpose[{texp, dexp}];
(* Above is my experimental data *) 

runlength = 2000; (* How many iterations for the simulation *)

(* Function to fit the data, this is called repeatedly in the \
   algorithm with different sets of parameters. *)
chi2[ka_Real, kb_Real, dexp_List, texp_List] :=
  {sol, A, B, c, t},
   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]]; 
Timing[chi2[0.1, 0.1, deep, texp]] (* How long does it take, what is the score *) 
(* Above is my function to fit the data, where A, B, c are the time dependent variables,
 and ka and kb are the model rate parameters that are adjusted to fit the data. *)  

(* Clearing stuff just to be safe *) 
 Clear[u, v] 

 (* 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}]; 
 (* End Storage arrays *) 

 range = 10.0;(* Limit initial guesses for parameter sets *) 
 weight = 1./10; 

 RandomSeed[AbsoluteTime[DateString[]]];(* Set Random Seed based on current time *) 

 ka[[1]] = RandomReal[{0.0, 200.}]; (* 1 random real number between 0.0 and 20.0 (or range  if specified) *) 
 kb[[1]] = RandomReal[{0.0, 20.}]; 
 least = 0.0004; (* Lower bound *) 
 max = 100.000;(* Upper bound *) 
 (* Give me an initial score for the first set of parameters above *) 
 decayscore[[1]] = chi2[ka[[1]], kb[[1]], dexp, texp]; 
 beta = 1/(2*0.05);(* This is a simulation temperature." *) 
 kar = RandomReal[{-1., 1.}, {runlength}]; 
 kbr = RandomReal[{-1., 1.}, {runlength}]; 
 u = RandomReal[{0., 1.}, {runlength}]; 

(******** Starting the Metropolis Algorithm ********) 

 tbl = 
   Table[(*Open Main For Loop*) 
    ka[[i]] = ka[[i - 1]] + weight*ka[[i - 1]]*kar[[i]]; (* Choosing new parameters *) 
    kb[[i]] = kb[[i - 1]] + weight*kb[[i - 1]]*kbr[[i]];(* Choosing new parameters *) 

  (*This is where the evaluation of the parameter sets begins*) 

  If[ (* Make sure the parameters are within a certain range *) 
   least < ka[[i]] && ka[[i]] < max &&  
   least < kb[[i]] && kb[[i]] < max, 

   decayscore[[i]] = chi2[ka[[i]], kb[[i]], dexp, texp]; 
   deltaE[[i]] = decayscore[[i]] - decayscore[[i - 1]]; 
   freeE[[i]] = Exp[-beta*deltaE[[i]]]; 

   v[[i]] = Min[1, freeE[[i]]]; 

   (******** Metropolis Selection Algorithm********) 

    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]]; 
    Sow[{True, decayscore[[i]], ka[[i]], kb[[i]], i}],(* Writing to bestParamSets *) 
    ka[[i]] = ka[[i - 1]]; 
    kb[[i]] = kb[[i - 1]]; 
   ka[[i]] = ka[[i - 1]]; 
   kb[[i]] = kb[[i - 1]] 
  (*Close Main For Loop//How long does this loop take to run?*) 
  {i, 2, runlength - 1}]][[2, 1]]; 
bestParamSets = Cases[tbl, {True, data__} -> {data}];) // AbsoluteTiming 
  • $\begingroup$ Are you aware that Mma starts counting with 1? With e.g. ka[[0]] = RandomReal[{0.0, 200.}]; you are replacing the head of ka (List) by a random number. This will undermine any code vectorization attempt. $\endgroup$ – Karsten 7. Nov 5 '14 at 11:15
  • $\begingroup$ It would be easier to help and prevent one form misinterpreting your code, if you would state your problem in a form like: this is my experimental data, this is the fit function with the following parameters and the variable .... Also, please explain why you can't use NonlinearModelFit. $\endgroup$ – Karsten 7. Nov 5 '14 at 11:28
  • $\begingroup$ @Karsten7. I have updated the code to begin counting at one, and have also edited the question per your suggestions. This simple example will converge via built in Mathematica functions. However the actual data and models that I am running are very complicated, and have proven to difficult to fit with the built in functions. Or cannot be programmed in a manner to satisfy the modeling conditions with the built in functions. $\endgroup$ – tarhawk Nov 5 '14 at 11:52
  • $\begingroup$ In this example there is a symbolic solution for A (35 E^(-ka t)). Are the actual models more complicated because there is no symbolic fit function you can use and you only have an InterpolatingFunction to compare your data to, or what makes them too difficult to fit with the build in functions? $\endgroup$ – Karsten 7. Nov 5 '14 at 12:17
  • $\begingroup$ @Karsten7. Yes, that is the case. The actual models can only be solved numerically, contain many different coupled odes, and have different variations on the initial conditions. There is also a large number of free parameters for each of the models. $\endgroup$ – tarhawk Nov 5 '14 at 12:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.