Help with NonLinearModelFit

Question

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

.

Answer 1

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

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

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

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

Tracking the iterations results in the following:

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.

Answer 2

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

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]

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