3
$\begingroup$

I have a question about fitting data with a model. I have a quite complicated model, based on Arctans, modified Bessel functions, etc, that I want to fit to my data, a table of couples (minimal example).

The first method that I used (brute force) with NonLinearModelFit works but return some warnings that I don't know how to avoid, and the fitting parameters, especially the last one are wrong.

Here is the code that I used:

ClearAll["Global`.*"];

g = 4.295676165008085`*^-6; 
y = a/b;

c[x_] := v/( (1 + x*y)*(1 + (x*y)^2) );
d0 = v*a^3;
e[x_] := d0 * (-ArcTan[x*y] + 0.5* Log[1 + (x*y)^2] + 
 Log[1 + (x*y)]);
f[x_] := Sqrt[g*e[x]/(x*a)]
(*h[x_] := Sqrt[3*g*h1/a]*(3*x)^2*(BesselI[0, 3*x]*BesselK[0, 3*x] - 
BesselI[1, 3*x]*BesselK[1, 3*x])*)
(*after reedit*)
h[x_] := 
Sqrt[3*g*h1/a]*(3*x)*Sqrt[(BesselI[0, 3*x]*BesselK[0, 3*x] - BesselI[1, 3*x]*BesselK[1, 3*x])]
j[x_] := (x*a*h[x]^2)/g 
(*k[x_] := Sqrt[g*k1/a]*x^2*(BesselI[0, x]*BesselK[0, x] - BesselI[1, x]*BesselK[1, x])*)
(*after reedit*)
k[x_]= 
Sqrt[g*k1/a]*x* Sqrt[(BesselI[0, x]*BesselK[0, x] - BesselI[1, x]*BesselK[1, x])]
l[x_] := (x*a*k[x]^2)/g 

(*free parameters of my fit : h1,k1,a,b,v*)

firstbin = 
  {{0.0745682, 5.77, 2.59}, {0.090855, 2.26, 3.66}, 
   {0.0815721, 13.1, 8.53}, {0.0505227, 7.89362, 1.16673}, 
   {0.0626294, 9.1, 10.4523}, {0.0772834, 3.89214, 5.4580}, 
   {0.084189, 3.7676, 0.714178}, {0.0927923, 1.7, 5.5}, 
   {0.0713416, 9.703, 2.34}, {0.0876698, 7.1589, 2.84105}, 
   {0.0376569, 8.85545, 5.53406}, {0.045388, 6.4866, 1.61059}, 
   {0.0559153, 4.67013, 1.75362}, {0.0801532, 9.81, 1.6915}, 
   {0.061584, 2.875, 1.90919}, {0.0930765, 4.5325456, 2.74357}, 
   {0.0832263, 7.306, 2.11707}, {0.0394981, 1.95511, 1.67028}, 
   {0.126419, 13.54, 3.29974}, {0.185348, 12.6713, 2.9279}, 
   {0.184662, 8.1504, 2.96985}, {0.154398, 8.72227, 3.84503}, 
   {0.199546,10.1779, 1.97452}};

z = Table[Log[10, firstbin[[i, 2]]], {i, 1, Length[firstbin]}];
rnorm1 = Table[firstbin[[i, 1]], {i, 1, Length[firstbin]}];
data1 = Table[{rnorm1[[i]], z[[i]]}, {i, 1, Length[z]}];

u[x_] = h[x] + f[x] + k[x];
fitLine1 = 
  NonlinearModelFit[
    data1, 
    {u[s], (k1/a > 0 && b > 0 && (a/b) > 0 && (h1/a) > 0)}, 
    {h1, k1, a, b, v}, 
    s]
fitLine1[s]

NonlinearModelFit::eit: The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 500 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {0.000037871,0.000885062,8.34782*10^-8}, is returned. >>

fitLine1["ParameterTable"]

FittedModel::constr: The property values {ParameterTable} assume an unconstrained model. The results for these properties may not be valid, particularly if the fitted parameters are near a constraint boundary. >>

This method works in 14.802 sec (determined with Timing), so it is quite long (I mean I think there is faster methods).

I tried with another method, but it gave me:

NonlinearModelFit::nrnum: The function value [...]is not a real number at {h1,k1,a,b,v} = {1.,1.,1.,1.,1.}. >>

Here is the code :

u[h1_?NumberQ, k1_?NumberQ, a_?NumberQ, b_?NumberQ, v_?NumberQ] :=  
  (u[h1, k1, a, b, v] = h[s] + f[s] + k[s])

fitLine1 = 
  NonlinearModelFit[
    data1, 
    {u[h1, k1, a, b, v][s], (k1/a > 0 && b > 0 && (a/b) > 0 && (h1/a) > 0)}, 
    {h1, k1, a, b, v}, 
    s]

I don't know why this doesn't work ( as I precised some conditions on the fit parameters), and if there is a way to make an optimal fit of my data?

PS: in the attached file, I put the table of parameters concerning my first attempt

fittable.

$\endgroup$
  • $\begingroup$ You have 6 data points and 6 parameters (5 fixed effects and the variance). Did you only use 6 data points as part of a minimal example (which needs to be less minimal) or do you really only have 6 data points? In addition, you see that all of the parameters are very near 1.0 which is what (I believe) Mathematica uses as default starting values. That's not a good sign. $\endgroup$ – JimB Oct 17 '15 at 19:16
  • $\begingroup$ well effectively, I agree that maybe the example is too minimal. I have actually more points. I will make a modeification in my question to include more points. The parameters are near 1, and this I don't understand actually... $\endgroup$ – Jean-Philippe Fontaine Oct 17 '15 at 21:22
  • 1
    $\begingroup$ The fact that all parameters are near 1 is just consequence of using the default start value (1) for all parameters by not specifying other start values. $\endgroup$ – Karsten 7. Oct 17 '15 at 21:54
  • $\begingroup$ I reduced each data point to a set of three quantities, where the last one is the uncertainty on the second one. $\endgroup$ – Jean-Philippe Fontaine Oct 17 '15 at 21:54
7
$\begingroup$

Update due to the formulas being updated in the question.

Using the updated equations we have the following:

(* Define equations *)
g = 4.295676165008085`*^-6;
d0[a_, v_] := v*a^3;
e[x_, a_, b_, v_] := 
  d0[a, v]*(-ArcTan[x*a/b] + (1/2)*Log[1 + (x*a/b)^2] + Log[1 + (x*a/b)]);
f[x_, a_, b_, v_] := Sqrt[g*e[x, a, b, v]/(x*a)]
h[x_, a_, h1_] := 
 Sqrt[3*g*h1/a]*(3*x)^2*(BesselI[0, 3*x]*BesselK[0, 3*x] - BesselI[1, 3*x]*BesselK[1, 3*x])
k[x_, a_, k1_] := 
 Sqrt[g*k1/a]*x^2*(BesselI[0, x]*BesselK[0, x] - BesselI[1, x]*BesselK[1, x])
u[x_, a_, b_, h1_, k1_, v_] = h[x, a, h1] + f[x, a, b, v] + k[x, a, k1];

When the function u is evaluated we have

Equation

We see that there are really only 4 parameters (rather than 5): h1/a, k1/a, v*a^2, and a/b. So the model is overparameterized and it should be fit with just 4 parameters.

But we're not quite done with highly linearly related pieces of the model. The two terms involving the Bessel functions are very much linearly related in the range of the data (x=0.04 to x=0.2). If predictions are what count rather than following some theoretical formula, then we can just remove one of the terms with the Bessel functions. And further, the Bessel function terms are both highly linearly related to the remaining term unless 10 < a/b < 60. Because all three of the pieces summed to make the final function are so highly linearly related, any fitting algorithm is going to have trouble finding a unique solution.

Below is a some code that can show the extent of the linear relationships:

(* Show the level of linear association with the 3 constituent parts of u2 *)
Manipulate[p4 = Exp[logp4]; GraphicsRow[{
   ParametricPlot[{x Sqrt[
      BesselI[0, x] BesselK[0, x] - BesselI[1, x] BesselK[1, x]], 
     Sqrt[(-ArcTan[p4 x] + Log[1 + p4 x] + 1/2 Log[1 + p4^2 x^2])/
     x]}, {x, 0.04, 0.2}, AspectRatio -> 1,
    Frame -> True, 
    FrameLabel -> {{"\!\(\*SqrtBox[FractionBox[\((\(-ArcTan[p4\\\ 
x]\) + Log[1 + p4\\\ x] + \*FractionBox[\(1\), \(2\)]\\\ Log[1 + 
\*SuperscriptBox[\(p4\), \(2\)]\\\ \*SuperscriptBox[\(x\), 
\(2\)]])\), \(x\)]]\)", 
       ""}, {"x \!\(\*SqrtBox[\(BesselI[0, x]\\\ BesselK[0, x] - 
BesselI[1, x]\\\ BesselK[1, x]\)]\)", 
       Dynamic[Style["p4 = " <> ToString[Exp[logp4]], Bold, Black, 
         Large]]}}, ImageSize -> 450],
   ParametricPlot[{x Sqrt[
      BesselI[0, x] BesselK[0, x] - BesselI[1, x] BesselK[1, x]], 
     x Sqrt[BesselI[0, 3 x] BesselK[0, 3 x] - 
       BesselI[1, 3 x] BesselK[1, 3 x]]}, {x, 0.04, 0.2},
    Frame -> True, 
    FrameLabel -> {{"x \!\(\*SqrtBox[\(BesselI[0, x]\\\ BesselK[0, x] 
- BesselI[1, x]\\\ BesselK[1, x]\)]\)", 
       ""}, {"x \!\(\*SqrtBox[\(BesselI[0, 3\\\ x]\\\ BesselK[0, 3\\\ 
x] - BesselI[1, 3\\\ x]\\\ BesselK[1, 3\\\ x]\)]\)", ""}}, 
    AspectRatio -> 1, ImageSize -> 400]}],
 {{logp4, Log[0.001], "p4 = a/b"}, Log[0.001], Log[10^8]}]

Linear relationships

In short, the pieces to the sum are all highly linearly related to each other and attempting to find the best linear sum of these pieces is somewhat fruitless (in Mathematica, R, SAS, or another other software package). Just fitting a simple function of the x values (maybe a low order polynomial?) would allow for just as good (if not better) predictions. But I would certainly like to be wrong in this case.

End of update.

You need to follow all of the advice of @JackLaVigne (especially the issue of the wildly different scales for the parameters). I know that you just wanted to give enough data to get the code working, however, it is likely that the model with 5 parameters is just too complex for 23 data points. Things might clear up with more data.

In addition, the ending coefficient estimates are found to be highly correlated which also tends to complicate convergence and makes interpreting coefficients individually a bit less certain. It might be that a re-parameterization is in order. For instance b, h1, and k1 are scaled by a. Maybe b/a, h1/a, and k1/a might be replaced with b, h1, and k1.

I modified the original code somewhat (and hopefully didn't change things incorrectly) and added in the use of the StepMonitor option to show the progress. I ran the model with 5,000 iterations and as you'll see below in the figures, convergence has yet to be reached with 5,000 iterations. The good news is that the convergence while slow at least appears to be steady: no wild fluctuations in the parameter estimates or of the mean square error.

g = 4.295676165008085`*^-6;

(* Define equations *)
d0[a_, v_] := v*a^3;
e[x_, a_, b_, v_] := 
  d0[a, v]*(-ArcTan[x*a/b] + 0.5*Log[1 + (x*a/b)^2] + Log[1 + (x*a/b)]);
f[x_, a_, b_, v_] := Sqrt[g*e[x, a, b, v]/(x*a)]
h[x_, a_, h1_] := 
 Sqrt[3*g*h1/a]*(3*x)^2*(BesselI[0, 3*x]*BesselK[0, 3*x] - BesselI[1, 3*x]*BesselK[1, 3*x])
k[x_, a_, k1_] := Sqrt[g*k1/a]*x^2*(BesselI[0, x]*BesselK[0, x] - BesselI[1, x]*BesselK[1, x])
u[x_, a_, b_, h1_, k1_, v_] = h[x, a, 10000000 h1] + f[x, a, b, v] + k[x, a, 10000000 k1];

(* Get the data in shape *)
firstbin = {{0.0745682, 5.77, 2.59}, {0.090855, 2.26, 
    3.66}, {0.0815721, 13.1, 8.53}, {0.0505227, 7.89362, 
    1.16673}, {0.0626294, 9.1, 10.4523}, {0.0772834, 3.89214, 
    5.4580}, {0.084189, 3.7676, 0.714178}, {0.0927923, 1.7, 
    5.5}, {0.0713416, 9.703, 2.34}, {0.0876698, 7.1589, 
    2.84105}, {0.0376569, 8.85545, 5.53406}, {0.045388, 6.4866, 
    1.61059}, {0.0559153, 4.67013, 1.75362}, {0.0801532, 9.81, 
    1.6915}, {0.061584, 2.875, 1.90919}, {0.0930765, 4.5325456, 
    2.74357}, {0.0832263, 7.306, 2.11707}, {0.0394981, 1.95511, 
    1.67028}, {0.126419, 13.54, 3.29974}, {0.185348, 12.6713, 
    2.9279}, {0.184662, 8.1504, 2.96985}, {0.154398, 8.72227, 
    3.84503}, {0.199546, 10.1779, 1.97452}};
z = Table[Log[10, firstbin[[i, 2]]], {i, 1, Length[firstbin]}];
rnorm1 = Table[firstbin[[i, 1]], {i, 1, Length[firstbin]}];
data1 = Table[{rnorm1[[i]], z[[i]]}, {i, 1, Length[z]}];
n = Length[data1[[All, 1]]]; (* Number of data points *)

(* Perform the fit and show the parameter values and resulting mean
   square error for each iteration *)
fitLine1 = 
 NonlinearModelFit[
  data1, {u[s, a, b, h1, k1, v], (a > 0 && b > 0 && k1 > 0 && h1 > 0 && v > 0)},
  {{h1, 2.65}, {k1, 2.2}, {a, 23.7}, {b, 0.005}, {v, 0.7}}, s, 
  MaxIterations -> 5000,  
  StepMonitor :> Print[{a, b, h1, k1, v, 
     Sum[(data1[[i, 2]] - u[data1[[i, 1]], a, b, h1, k1, v])^2/(n - 5), {i, n}]}]];

(* Show some of the fit statistics *)
fitLine1["BestFitParameters"]
(* {h1 -> 2.41488, k1 -> 3.74096, a -> 24.2358, b -> 0.0076914, v -> 0.755862} *)
fitLine1["EstimatedVariance"]
(* 0.073205 *)
fitLine1["CorrelationMatrix"] // MatrixForm

Correlation matrix

(* Plot the data and the fitted curve for the last iteration *)
Show[ListPlot[data1], Plot[fitLine1[s], {s, 0.04, 0.2}], Frame -> True]

Model fit

Tracking the iterations results in the following:

Iteration progress

It appears there is no leveling off and there's a long way to go. But, again, this might clear up with a larger number of sample points. And once you do get convergence, I'd look at the residuals to see if modifications in the formula and/or error structure might be warranted as it looks like there's more variability for low values of s than high values of s.

$\endgroup$
  • $\begingroup$ +1 for scaling the function so that the parameters are in the same ball park. The b parameter should probably also be scaled. I very much like defining the functions to include all the parameters (e.g., u[x_, a_, b_, h1_, k1_, v_]. I have validated your modification of the original code find it to be correct. $\endgroup$ – Jack LaVigne Oct 18 '15 at 10:43
  • $\begingroup$ Thank you very much for your two complete replies. I didn't react for the moment in order to understand the steps in the code you proposed me. I saw that I made 2 mistakes in my fitting, concerning the definition of h and k, I will edit my post. However the fit isn't superimposed to my datas in Manipulate. Is it because of the starting values of the parameters that I have to change manually until getting something? is there a way to do that, or I just choose random (but physically meaningful) values? $\endgroup$ – Jean-Philippe Fontaine Oct 18 '15 at 14:07
  • $\begingroup$ @Jean-PhilippeFontaine Yes, you will have to adjust parameters accordingly. Use 10^5 rather than 10^7 for the multiplier of h1 and k1. $\endgroup$ – Jack LaVigne Oct 18 '15 at 14:30
  • $\begingroup$ @ Jim Baldwin, I looked at the correlation matrix for my model (the reparametrized one that you proposed, and the other one), and I found coefficients all close to 1 (a bit less for the reparametrized model, between 0.91 and 0.99 instead of 0.99). So it seems that the fit parameters, even after reparametrization are still correlated. $\endgroup$ – Jean-Philippe Fontaine Oct 18 '15 at 23:54
6
$\begingroup$

I studied this problem and have some comments.

Equations and data from your problem

g = 4.295676165008085`*^-6;
y = a/b;
c[x_] := v/((1 + x*y)*(1 + (x*y)^2))
d0 = v*a^3;
e[x_] := d0*(-ArcTan[x*y] + 0.5*Log[1 + (x*y)^2] + Log[1 + (x*y)])
f[x_] := Sqrt[g*e[x]/(x*a)]
h[x_] := Sqrt[3*g*h1/a]*(3*x)^2*(BesselI[0, 3*x]*BesselK[0, 3*x] - 
    BesselI[1, 3*x]*BesselK[1, 3*x])
j[x_] := (x*a*h[x]^2)/g
k[x_] := Sqrt[g*k1/a]*
  x^2*(BesselI[0, x]*BesselK[0, x] - BesselI[1, x]*BesselK[1, x])
l[x_] := (x*a*k[x]^2)/g
u[x_] = h[x] + f[x] + k[x]

firstbin = {{0.0745682, 5.77, 2.59}, {0.090855, 2.26, 
    3.66}, {0.0815721, 13.1, 8.53}, {0.0505227, 7.89362, 
    1.16673}, {0.0626294, 9.1, 10.4523}, {0.0772834, 3.89214, 
    5.4580}, {0.084189, 3.7676, 0.714178}, {0.0927923, 1.7, 
    5.5}, {0.0713416, 9.703, 2.34}, {0.0876698, 7.1589, 
    2.84105}, {0.0376569, 8.85545, 5.53406}, {0.045388, 6.4866, 
    1.61059}, {0.0559153, 4.67013, 1.75362}, {0.0801532, 9.81, 
    1.6915}, {0.061584, 2.875, 1.90919}, {0.0930765, 4.5325456, 
    2.74357}, {0.0832263, 7.306, 2.11707}, {0.0394981, 1.95511, 
    1.67028}, {0.126419, 13.54, 3.29974}, {0.185348, 12.6713, 
    2.9279}, {0.184662, 8.1504, 2.96985}, {0.154398, 8.72227, 
    3.84503}, {0.199546, 10.1779, 1.97452}};

The data was formed via:

z = Log[10, firstbin[[All, 2]]];
rnorm1 = firstbin[[All, 1]];
data1 = Transpose[{rnorm1, z}];

You define c[x], l[x] and j[x] but they are not used anywhere. This is suspicious.

Data and model

I defined a local function that took the parameters as arguments. The expression was derived by evaluating u[s].

func[s_, h1_, k1_, a_, b_, v_] := 
 0.0020726 Sqrt[k1/a]
    s^2 (BesselI[0, s] BesselK[0, s] - BesselI[1, s] BesselK[1, s]) + 
  0.032308656 Sqrt[h1/a] s^2
    (BesselI[0, 3 s] BesselK[0, 3 s] - 
     BesselI[1, 3 s] BesselK[1, 3 s]) + 
  0.0020726 Sqrt[(
   a^2 v (-ArcTan[(a s)/b] + Log[1 + (a s)/b] + 
      0.5` Log[1 + (a^2 s^2)/b^2]))/s]

There are three terms: the first has k1, the second h1 and the third v and b. The parameter a is in all of the terms.

I used manually adjusted the parameters until I could get it to plot on the same scale as the data using Manipulate. Note that I had to use some huge numbers as well as very small numbers. Edit: Used parameters derived by Jim Baldwin

Note also that it doesn't match the data very well.

Manipulate[
 Grid[{
   {"k[s]", 
    0.002072601303919325` Sqrt[k1/a]
      s^2 (BesselI[0, s] BesselK[0, s] - BesselI[1, s] BesselK[1, s])},
   {"h[s]", 
    0.032308656859995975` Sqrt[h1/a]
      s^2 (BesselI[0, 3 s] BesselK[0, 3 s] - 
       BesselI[1, 3 s] BesselK[1, 3 s])},
   {"f[s]", 
    0.002072601303919325` \[Sqrt](1/
        s a^2 v (-ArcTan[(a s)/b] + Log[1 + (a s)/b] + 
          0.5` Log[1 + (a^2 s^2)/b^2]))},

   {"u[s]", func[s, h1, k1, a, b, v]},
   {Show[
     Plot[func[x, h1, k1, a, b, v],
      {x, 0.01, 0.2},
      PlotStyle -> Black,
      PlotRange -> {{0.035, 0.202}, {-0.01, 1.07}},
      AxesOrigin -> {0.035, -0.01}
      ],
     ListPlot[data1,
      PlotStyle -> Red
      ],
     ImageSize -> 300
     ], SpanFromLeft}
   }],
 {{s, 0.18}, 0, 0.2, Appearance -> "Open"},
 {{h1, 2.4*10^7}, 1.*10^7, 9.*10^7., Appearance -> "Open"},
 {{k1, 3.7*10^7}, 1*10^7, 9*10^7., Appearance -> "Open"},
 {{a, 24.}, 10, 30, Appearance -> "Open"},
 {{b, 0.0077}, 0, 0.01, Appearance -> "Open"},
 {{v, 0.76}, 0, 1, Appearance -> "Open"}
 ]

Mathematica graphics

A Way Forward

A good way to test optimization problems is to create synthetic data where you know the values of the parameters and make sure you have the syntax correct and that the spread in the data is sufficient to determine the parameters.

Some times one adds noise to the synthetic data. In this example I will not add noise.

fakeData = Table[{s, func[s, 2., 5., 0.1, 0.1, 10.]}, {s, 0.01, 0.20, 0.01}]

ListPlot[fakeData]

Mathematica graphics

I only had to make small adjustments to your non linear model expression. You had some constraints involving ratios. If you study those constraints you will that they are equivalent to defining the parameters to be positive.

You had the constraints surrounded by parenthesis. I believe they need to be enclosed by curly brackets.

Below is the problem using the fake data.

fitLine1 =  NonlinearModelFit[fakeData, 
  {u[s], {h1 > 0 && k1 > 0 && a > 0 && b > 0}},
  {h1, k1, a, b, v}, s]

It runs quickly but produces a poor result.

fitLine1["BestFitParameters"]

{h1 -> 0.988858, k1 -> 0.998287, a -> 1.27338, b -> 0.866405, v -> 16.2319}

Now try it again setting good starting values.

fitLine1 =  NonlinearModelFit[fakeData, 
  {u[s], {h1 > 0 && k1 > 0 && a > 0 && b > 0}},
  {{h1, 2}, {k1, 5.}, {a, 0.1}, {b, 0.1}, {v, 10}}, s]

fitLine1["BestFitParameters"]

{h1 -> 2.00184, k1 -> 5.00179, a -> 0.101587, b -> 0.101719, v -> 14.9816}

Most of the parameters match the input with the exception of v. If we have problems reproducing the input parameters it tells us that we certainly will have problems with real data.

It also shows that we must have reasonable starting values for any hope at all.

Simplifying the model

Jim Baldwin showed that there was terrible cross-correlation between the parameters. The optimizer is spinning its wheels because we have given it too many degrees of freedom.

A new function is defined were we have removed the parameter a and redefined the remaining parameters as follows:

Sqrt[h1/a] -> H1, Sqrt[k1/a] -> K1 ,  Sqrt[(a^2) v ] -> V, a/b -> B

Note that the new parameters have the square root applied and that B is proportional to the inverse of b.

The new function is:

U[s_, B_, H1_, K1_, V_] = 
 0.0020726 K1 s^2 (BesselI[0, s] BesselK[0, s] - 
     BesselI[1, s] BesselK[1, s]) + 
  0.0323087 H1 s^2 (BesselI[0, 3 s] BesselK[0, 3 s] - 
     BesselI[1, 3 s] BesselK[1, 3 s]) + 
  0.002072 V Sqrt[(-ArcTan[B s] + Log[1 + B s] + 
     0.5` Log[1 + B^2 s^2])/s]

Now we fit it. The starting values need to be modified to be consistent with the new parameters.

{fitLine, steps} = 
  Reap[NonlinearModelFit[
    data1, {U[s, B, H1, K1, 
      V], {B > 0 && K1 > 0 && H1 > 0 && V > 0}}, {{H1, 1050.}, {K1, 
      960.}, {B, 4740.}, {V, 20.}}, s, 
    StepMonitor :> 
     Sow[{H1, K1, B, V, 
       Sum[(data1[[i, 2]] - U[data1[[i, 1]], B, H1, K1, V])^2/(n - 
           5), {i, n}]}]]
   ]

This now converges nicely after about 60 iterations.

fitLine2["BestFitParameters"]

{H1 -> 737.279, K1 -> 2676.99, B -> 3579.65, V -> 21.4477}

We need to convert the parameters back to the original space. This gives

{h1 -> 5.436*10^5, k1 -> 7.166*10^6, a -> 1., b -> 0.0002794, v -> 460.}

Now we plot the data vs model

Mathematica graphics

$\endgroup$
  • $\begingroup$ the problem when we go back to the original space is that the system of equations is surdimensionned, so we have to set one parameter to 1 (a in your case). Also I am not that familiar with Reap Sow and I don't know what is the expression you calculated in Sow, is it the variance ? $\endgroup$ – Jean-Philippe Fontaine Oct 18 '15 at 16:25
  • $\begingroup$ oh sorry I didn't pay attention carefully, it is defined in the reply of Jim Baldwin $\endgroup$ – Jean-Philippe Fontaine Oct 18 '15 at 18:26
  • $\begingroup$ @ Jack LaVigne; I tried to fit with NonLinearModelFit (wrapped into a Reap..Sow) but I still got a message "NonlinearModelFit::eit: "The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 500 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {2.62704*10^-9,0.000441404,9.19079*10^-10}, is returned. " Also as I mentionned above, doing the reverse transformation to get the values of the initial paremeters, we need to fix one to get the others ; otherwise the system of equations is surdimensionned. $\endgroup$ – Jean-Philippe Fontaine Oct 18 '15 at 21:56
  • $\begingroup$ I am using datas that I provided in my question. Actually it gives answer but with this warning. $\endgroup$ – Jean-Philippe Fontaine Oct 18 '15 at 23:14
  • $\begingroup$ @Jean-PhlippeFontaine You can skip the fakeData portion but you should be able to copy and paste from the start with the modified h and k functions. You will need to redo the starting points for the simplified model. I used {{H1, 160.}, {K1, 200.}, {B, 1.}, {V, 3240.}}. It converged quickly. Good luck. $\endgroup$ – Jack LaVigne Oct 19 '15 at 0:32

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.