NonlinearModelFit returning imaginary numbers after fitting real data

Question

I am getting problems when I try to use NonlinearModelFit to fit a squareroot-like function to a set of real data points.

Here is the function:

\begin{equation} \left(-\frac{2 \sum _{i=0}^n c_i q^i}{\sqrt{\left(\sum _{i=0}^M b_i q^i\right){}^2-4 \left(a_1 q+1\right) \sum _{i=0}^n c_i q^i}+\sum _{i=0}^M b_i q^i}\right)^2 \end{equation}

where $M,\,n\in \mathbb{N}$.

I can't put constraints on NonLinearModelFit because I need the errors, and I can't try to modify the function with Re or Conjugate. Here is the function for M=2 and n=2:

Model[q_]:=(4*(Subscript[c, 0] + q*Subscript[c, 1] + q^2*Subscript[c, 2])^2)/(Subscript[b, 0] + q*Subscript[b, 1] + 2*Subscript[b,2]+Sqrt[(Subscript[b, 0] + q*Subscript[b, 1] + q^2*Subscript[b, 2])^2 - 4*(1 + q*Subscript[a, 1])*(Subscript[c, 0] + q*Subscript[c, 1] +  q^2*Subscript[c, 2])])^2

Here is the NonLinearModelFit:

NonlinearModelFit[Points,Model[q], {Subscript[a, 1], Subscript[b, 0], Subscript[b, 1], Subscript[b, 2], Subscript[c, 0], Subscript[c, 1], Subscript[c, 2]}, q]

Here is the data:

{{{0.0491483, 0.85451}, {0.00306566, 0.989755}, {0.156523, 0.634745}, {0.19322, 0.579874}, {0.107026, 0.723105}, {0.159402, 0.630146}, {0.126612, 0.685905}, {0.041859, 0.873807}, {0.202128, 0.567694}, {0.00429277, 0.985704}, {0.0438898, 0.868359}, {0.202524, 0.567161}, {0.164643, 0.621908}, {0.0176414, 0.943355}, {0.0887417, 0.760861}, {0.247813, 0.511176}, {0.2382, 0.522303}, {0.215778, 0.549804}, {0.140331, 0.661648}, {0.0450278, 0.865331}, {0.079667, 0.780801}, {0.200001, 0.570565}, {0.13433, 0.672086}, {0.130211, 0.679405}, {0.244809, 0.514612}, {0.181252, 0.596908}, {0.104507, 0.728124}, {0.171583, 0.611263}, {0.0359265, 0.890045}, {0.187224, 0.58831}, {0.225738, 0.53731}, {0.0723277, 0.797555}, {0.197513, 0.573953}, {0.10055, 0.736123}, {0.141649, 0.65939}, {0.0926702, 0.752483}, {0.0598971, 0.827301}, {0.176524, 0.60386}, {0.0342446, 0.894738}, {0.110842, 0.715605}, {0.0116676, 0.961926}, {0.163328, 0.623958}, {0.243646, 0.515952}, {0.0561636, 0.836589}, {0.203855, 0.565379}, {0.0780196, 0.784511}, {0.196865, 0.574841}, {0.175081, 0.606007}, {0.0955304, 0.746478}, {0.0320406, 0.900951}, {1.5931, 0.0890937}, {1.0712, 0.146087}, {1.75196, 0.078463}, {2.65825, 0.0434119}, {0.827541, 0.194791}, {1.41378, 0.104037}, {2.60719, 0.0446751}, {1.56054, 0.0915489}, {0.647142, 0.249531}, {2.86556, 0.0388113}, {2.65311, 0.0435367}, {2.67742, 0.0429519}, {0.834751, 0.193004}, {1.00899, 0.156551}, {1.5725, 0.0906352}, {2.34786, 0.0520663}, {0.270901, 0.48594}, {0.871762, 0.184223}, {1.24283, 0.122258}, {1.30046, 0.115596}, {0.489775, 0.319958}, {0.501975, 0.31344}, {1.51357, 0.0952797}, {2.58615, 0.0452123}, {2.08624, 0.0616378}, {1.91234, 0.0696033}, {0.646257, 0.249857}, {1.55824, 0.0917261}, {2.11203, 0.0605774}, {0.580136, 0.276214}, {4.01178, 0.0230932}, {10.7232, 0.00449478}, {11.2478, 0.00413793}, {8.94916, 0.00613476}, {6.91823, 0.00948793}, {5.38288, 0.0143832}, {8.16001, 0.00718018}}}

Remind, put constraints in the model increase a lot the error.

Answer 1

I copied your data and called it points. In order to plot it as a line I sorted it.

points = SortBy[points, First]

There have been many posts on this site indicating that a Subscript is not a symbol. It can wreak havoc when used in Mathematica functions that expect symbols and should be avoided.

Mathematica uses uppercase symbols (e.g., Pi, NonlinearModelFit) for system symbols. Users are encouraged to use lowercase symbols when defining functions.

Thus it will be easier to write your model as:

model[q_] := (4*(c0 + q*c1 + q^2*c2)^2)/(b0 + q*b1 + q^2*b2 + 
     Sqrt[(b0 + q*b1 + q^2*b2)^2 - 4*(1 + q*a1)*(c0 + q*c1 + q^2*c2)])^2

Whatever approach one uses with a problem as complicated as this, the optimization function is going to need a reasonable starting point.

Manipulate can be used to vary the parameters and select values that give a reasonable fit. This is not an easy step but may be (depending upon the problem) required in order to get a good fit.

Manipulate[
 Show[
  Plot[
   (4*(c0 + q*c1 + q^2*c2)^2)/(b0 + q*b1 + q^2*b2 + 
       Sqrt[(b0 + q*b1 + q^2*b2)^2 - 
         4*(1 + q*a1)*(c0 + q*c1 + q^2*c2)])^2, 
   {q, 0.003, 12},
   PlotRange -> {{0, 12}, {0, All}},
   PlotStyle -> Blue
   ],
  ListLinePlot[points, PlotStyle -> Red],
  ImageSize -> 400
  ],
 {{a1, 0.2, "a1"}, 0, 5, Appearance -> "Open"},
 {{b0, 2, "b0"}, 0, 5, Appearance -> "Open"},
 {{b1, 0.55, "b1"}, 0, 5, Appearance -> "Open"},
 {{b2, 2.4, "b2"}, 0, 5, Appearance -> "Open"},
 {{c0, 1, "c0"}, -5, 5, Appearance -> "Open"},
 {{c1, 0.3, "c1"}, -5, 5, Appearance -> "Open"},
 {{c2, 0.4, "c2"}, -5, 5, Appearance -> "Open"}
 ]

This produced the following plot

The parameters from above were used as the starting point for NonlinearModelFit

nlm = NonlinearModelFit[points,
  (4*(c0 + c1*q + c2*q^2)^2)/(b0 + b1*q + b2*q^2 + 
      Sqrt[(b0 + b1*q + b2*q^2)^2 - 4*(1 + a1*q)*(c0 + c1*q + c2*q^2)])^2,
  {
   {a1, 0.2},
   {b0, 2},
   {b1, 0.55},
   {b2, 2.4},
   {c0, 1},
   {c1, 0.3},
   {c2, 0.4}
   }, q]

The symbol nlm is known as a FittedModel. It can be used as a function and can produce many diagnostics. The parameters can be seen by using "BestFitParameters" as the argument to nlm.

nlm["BestFitParameters"]

{a1 -> 1.42754, b0 -> 3.99063, b1 -> 6.07, b2 -> 1.68324, 
 c0 -> 2.99061, c1 -> 1.28942, c2 -> 0.00640726}

The fit is excellent.

Show[
 ListPlot[points, PlotStyle -> Red],
 Plot[
  nlm[q], 
  {q, 0.003, 12},
  PlotRange -> {{0, 12}, All},
  PlotStyle -> Blue
  ]
 ]

A word of caution about the result. A pretty printed version of the "CorrelationMatrix" is displayed below.

There are many parameters that have a high degree of correlation. This means that the problem is very poorly resolved.