Updating procedural fitting algorithm to more efficient style?

Question

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!

Thanks,

Clear["Global`*"]
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] :=
 Module[
  {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[bestParamSets] 
 Clear[deltaE] 
 Clear[freeE] 
 Clear[ka] 
 Clear[kb] 
 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 = 
  Reap[ 
   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********) 
   (**************************************************) 

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