0
$\begingroup$

I am fitting some model function to five sets of data through the following short code 

Needs["ErrorBarPlots`"]

MasData1 = {{{89, 6.7}, ErrorBar[1.272]}, {{112, 7.9}, 
ErrorBar[1.221]}, {{141, 9.3}, ErrorBar[1.697]}}; 
MasData2 = {{{83.9, 4.04}, ErrorBar[0.7754]}, {{114.1, 5.29}, 
ErrorBar[1.086]}, {{144.2, 6.1}, ErrorBar[1.681]}}; 
MasData3 = {{{62, 16.6}, ErrorBar[2.6172]}, {{85, 20.7}, 
ErrorBar[3.0809]}, {{108, 21.9}, ErrorBar[3.0647]}, {{135, 25.8}, 
ErrorBar[3.9115]}, {{183, 33.2}, ErrorBar[5.993]}, {{83.9, 14.5}, 
ErrorBar[2.772]}, {{114.1, 24.7}, 
ErrorBar[4.5875]}, {{144.2, 24.1}, ErrorBar[6.5756]}}; 
MasData4 = {{{53.3, 25.1}, ErrorBar[3.5489]}, {{83.9, 30}, 
ErrorBar[4.2095]}, {{114.1, 41.5}, 
ErrorBar[6.1404]}, {{144.2, 45}, ErrorBar[9.6243]}, {{57, 27.4}, 
ErrorBar[3.6056]}, {{80, 36.7}, ErrorBar[7.9925]}, {{101, 43}, 
ErrorBar[6.6138]}, {{128, 48.8}, ErrorBar[9.1001]}, {{180, 61.1}, 
ErrorBar[10.5575]}}; 
MasData5 = {{{44.8, 47.5}, ErrorBar[4.0]}, {{54.8, 50.1}, 
ErrorBar[4.2]}, {{64.8, 61.7}, ErrorBar[5.1]}, {{74.8, 64.8}, 
ErrorBar[5.5]}, {{84.9, 75}, ErrorBar[6.2]}, {{94.9, 81.2}, 
ErrorBar[6.7]}, {{104.9, 85.3}, ErrorBar[7.1]}, {{119.5, 94.5}, 
ErrorBar[7.5]}, {{144.1, 101.5}, ErrorBar[8.3]}, {{144.9, 101.9}, 
ErrorBar[10.9]}, {{162.5, 117.8}, 
ErrorBar[12.8]}, {{177.3, 130.2}, 
ErrorBar[13.4]}, {{194.8, 147.7}, 
ErrorBar[17.1]}, {{219.6, 137.4}, 
ErrorBar[20.1]}, {{244.8, 176.6}, 
ErrorBar[20.3]}, {{267.2, 178.7}, 
ErrorBar[21.1]}, {{292.3, 200.4}, ErrorBar[29.1]}, {{60, 55.8}, 
ErrorBar[4.838]}, {{80, 66.6}, ErrorBar[7.280]}, {{100, 73.4}, 
ErrorBar[6.426]}, {{120, 86.7}, ErrorBar[7.245]}, {{140, 104}, 
ErrorBar[12.083]}, {{160, 110}, ErrorBar[16.279]}, {{42.5, 43.8}, 
ErrorBar[3.482]}, {{55, 57.2}, ErrorBar[3.980]}, {{65, 62.5}, 
ErrorBar[4.614]}, {{75, 68.9}, ErrorBar[5.197]}, {{85, 72.1}, 
ErrorBar[5.523]}, {{100, 81.9}, ErrorBar[5.368]}, {{117.5, 95.7}, 
ErrorBar[6.277]}, {{132.5, 103.9}, ErrorBar[6.912]}, {{155, 115}, 
ErrorBar[7.920]}, {{185, 129.1}, ErrorBar[9.192]}, {{215, 141.7}, 
ErrorBar[10.666]}, {{245, 140.3}, ErrorBar[14.526]}, {{275, 189}, 
ErrorBar[24.274]}, {{49, 39.2}, ErrorBar[10]}, {{86, 75.7}, 
ErrorBar[14.414]}, {{167, 118}, ErrorBar[22.828]}, {{43.2, 50.7}, 
ErrorBar[1.5]}, {{50, 59.5}, ErrorBar[1.4]}, {{57.3, 61.8}, 
ErrorBar[1.9]}, {{65.3, 67.6}, ErrorBar[1.7]}, {{73.9, 72.4}, 
ErrorBar[1.9]}, {{83.2, 79.9}, ErrorBar[2.3]}, {{93.3, 84.4}, 
ErrorBar[2.1]}, {{104.3, 86.7}, ErrorBar[2.7]}, {{47.9, 55.4}, 
ErrorBar[2.1]}, {{68.4, 66.4}, ErrorBar[2.9]}};

gamma = 5.55*^-6;
p = 3.1;
alphaem = 1/137;
alphas1 = 0.2478239650859146
alphas2 = 0.2390601794032581
alphas3 = 0.266809332867253
alphas4 = 0.2796636697708153
alphas5 = 0.295705
Rg = 2^(2*(a + b*Log[q/0.45]) + 3)/Sqrt[Pi]*Gamma[(a + b*Log[q/0.45]) + 5/2]/ Gamma[(a + b*Log[q/0.45]) + 4];

F1[w_] = 3.89379*^5*1/(4.9 + 4*0.06*Log[w/90])*gamma*
p^3*(Pi)^3/48/
 alphaem*(alphas1/(6.4025)^2*Rg*
   NN*((j + p^2)/(j + w^2))^(-a - 
      b*Log[6.4025/0.45]))^2*2.665*(1 + (Pi)^2/
    4*(a + b*Log[q/0.45])^2) /. {q -> 6.4025, j -> 16};
F2[w_] = 3.89379*^5*1/(4.9 + 4*0.06*Log[w/90])*gamma*
p^3*(Pi)^3/48/
 alphaem*(alphas2/(8.0025)^2*Rg*
   NN*((j + p^2)/(j + w^2))^(-a - 
      b*Log[8.0025/0.45]))^2*3.331*(1 + (Pi)^2/
    4*(a + b*Log[q/0.45])^2)  /. {q -> 8.0025, j -> 22.4};
F3[w_] = 3.89379*^5*1/(4.9 + 4*0.06*Log[w/90])*gamma*
p^3*(Pi)^3/48/
 alphaem*(alphas3/(4.1525)^2*Rg*
   NN*((j + p^2)/(j + w^2))^(-a - 
      b*Log[4.1525/0.45]))^2*1.728*(1 + (Pi)^2/
    4*(a + b*Log[q/0.45])^2) /. {q -> 4.1525, j -> 7};
F4[w_] = 3.89379*^5*1/(4.9 + 4*0.06*Log[w/90])*gamma*
p^3*(Pi)^3/48/
 alphaem*(alphas4/(3.2025)^2*Rg*
   NN*((j + p^2)/(j + w^2))^(-a - 
      b*Log[3.2025/0.45]))^2*1.323*(1 + (Pi)^2/
    4*(a + b*Log[q/0.45])^2)  /. {q -> 3.2025, j -> 3.1};
F5[w_] = 3.89379*^5*1/(4.9 + 4*0.06*Log[w/90])*gamma*
p^3*(Pi)^3/48/
 alphaem*(alphas5/(2.4025)^2*Rg*
   NN*((j + p^2)/(j + w^2))^(-a - 
      b*Log[2.4025/0.45]))^2*1*(1 + (Pi)^2/
    4*(a + b*Log[q/0.45])^2)  /. {q -> 2.4025, j -> 0};

 chisq1 = Sum[((MasData1[[k]][[1]][[2]] - F1[MasData1[[k]][[1]][[1]]])/
  MasData1[[k]][[2]][[1]])^2, {k, 1, Length[MasData1]}];
 chisq2 = Sum[((MasData2[[k]][[1]][[2]] - F2[MasData2[[k]][[1]][[1]]])/
  MasData2[[k]][[2]][[1]])^2, {k, 1, Length[MasData2]}];
 chisq3 = Sum[((MasData3[[k]][[1]][[2]] - F3[MasData3[[k]][[1]][[1]]])/
  MasData3[[k]][[2]][[1]])^2, {k, 1, Length[MasData3]}];
 chisq4 = Sum[((MasData4[[k]][[1]][[2]] - F4[MasData4[[k]][[1]][[1]]])/
  MasData4[[k]][[2]][[1]])^2, {k, 1, Length[MasData4]}];
 chisq5 = Sum[((MasData5[[k]][[1]][[2]] - F5[MasData5[[k]][[1]][[1]]])/
  MasData5[[k]][[2]][[1]])^2, {k, 1, Length[MasData5]}];

  chisq = chisq1 + chisq2 + chisq3 + chisq4 + chisq5;
  rr = FindMinimum[chisq, a, b, NN]

  Show[Plot[{(F1[w] /. {a -> 0.05, b -> 0.079, NN -> 1.21}), (F2[
  w] /. {a -> 0.05, b -> 0.079, NN -> 1.21}), (F3[
  w] /. {a -> 0.05, b -> 0.079, NN -> 1.21}), (F4[
  w] /. {a -> 0.05, b -> 0.079, NN -> 1.21}), (F5[
  w] /. {a -> 0.05, b -> 0.079, NN -> 1.21})}, {w, 10, 300},Frame -> True, PlotRange -> {0, 210}, BaseStyle -> FontSize -> 18, PlotStyle -> {{Thick, Blue}, {Thick, Green}, {Thick, Red}, {Thick, Orange}, {Thick, Purple}}, PlotLegends -> Placed[{"DataSet 1", "DataSet 2", "DataSet 3", "DataSet 4", "DataSet 5"}, Scaled[{0.1, 0.75}]], Axes -> False],ErrorListPlot[{MasData1, MasData2, MasData3, MasData4, MasData5},PlotStyle -> {{Thick, Blue}, {Thick, Green}, {Thick, Red}, {Thick,Orange}, {Thick, Purple}}]]

Basically the whole procedure via a chi square minimisation attempts to extract the best fit parameters for the coefficients $NN,a,b$ within the model function. However, if one runs the code, the chi square returned is around a ridiculous 16000 not including division per degree of freedom. I am repeating some exercise so I know that I should find that the best fit parameters are $NN=1.21, a=0.05$ and $b=0.079$. I plot the model functions $Fi$ with these constants hardwired and as can be seen from the final plot, the fit is very well. So, why won't mathematica give me these constants? Is there a fault in my code? Thanks in advance! This error has evaded me for one week now.

|improve this question
$\endgroup$
1
$\begingroup$

FindMinimum looks for a local minima. Minimize attempts to find a global minimum.

rr = Minimize[chisq, {a, b, NN}] yields

{57.767, {a -> 0.0407385, b -> 0.076216, NN -> 1.38399}}

with plot fit:

fitplot

This indeed seems smaller than the $\chi^2$ of the older values you were aiming for: chisq /. {NN -> 1.21, a -> 0.05, b -> 0.079} $\mapsto~ 114.679$

Why did FindMinimum fail?

It started out in the weeds, and did manage to find somewhat of local valley, at first glance (See below) there's no nice beautiful descent to the global minima (although maybe if you squint you can see a tilt?). If you look at $\chi^2$ in the small region around where it exited out, you'll see this:

Plot3D[chisq /. a -> 0.7177534626280686`, {b, .8*.41095032157812883`, 
  1.20*.41095032157812883`}, {NN, 3.4843459401817767`*^-6 *.8, 
  3.4843459401817767`*^-6*1.2}, PlotStyle -> Directive[Opacity[.8]], 
 ClippingStyle -> {None}, PlotRange -> {15000, 20000}]

badValley

Compare with the nice smooth rolling valley near the global minima:

Plot3D[chisq /. a -> .04 , {b, .07, .09}, {NN, 1.1, 1.6}, 
 PlotRange -> {50, 150}, PlotPoints -> 50, 
 PlotStyle -> Directive[Opacity[.8]], ClippingStyle -> {None}]

goodValley

FindMinimum can be triumphant

To see why it pooped out (it didn't seem to be swamped with precision glitches), I decided to look at the rate of convergence:

ListPlot[Last[
  Reap[FindMinimum[chisq, a, b, NN, EvaluationMonitor :> Sow[chisq], 
    MaxIterations -> 100]]], PlotRange -> {15000, 20000}]

Plotting $\chi^2$ vs steps:

chiSqvsSteps

It looks like FindMinimum was simply frustrated that it was taking so long, but it didn't look like it was done minimizing. What happens when we allow some more time?

ListPlot[Last[
  Reap[FindMinimum[chisq, a, b, NN, EvaluationMonitor :> Sow[chisq], 
    MaxIterations -> 1000]]], PlotRange -> {0, 20000}]

Low and behold:

enter image description here

we do roll exactly to:

{57.767, {a -> 0.0407385, b -> 0.076216, NN -> 1.38399}}

So in this case default method of Find-Minimum was able to find the global minimum when given enough time!

|improve this answer
$\endgroup$
  • $\begingroup$ Oh! Thanks very much! I just wonder in what sense is this 16000 chi square local minimum really a local minimum because as you probably saw from the graph that could be produced in my OP, the fit is not satisfactory at all. Thanks! $\endgroup$ – CAF Aug 15 '17 at 21:20
  • $\begingroup$ FindMinimum probably failed to converge -- if it started out in the wilderness I expect it got nowhere in 100 iterations and just quit out. Indeed trying it now I get a FindMinimum::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. response. $\endgroup$ – John Joseph M. Carrasco Aug 15 '17 at 21:23
  • $\begingroup$ I see. So does iteration mean like a scan of the multidimensional $\chi^2, NN,a,b$ space? So the 16000 just like represents a random point in the $\left\{\chi^2, NN,a,b\right\}$ space or something? $\endgroup$ – CAF Aug 15 '17 at 21:28
  • $\begingroup$ Yes to both. I'll add some notes to the answer w/ some plots. $\endgroup$ – John Joseph M. Carrasco Aug 15 '17 at 21:40
  • $\begingroup$ I see! if FindMinimum's job is to find a local minimum and 16000 looks like a local minimum to me from your plots then what does it mean that it failed to converge? $\endgroup$ – CAF Aug 15 '17 at 21:57

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.