# NonLinearFit with “complicated” integral model

I have a set of data $y_i (z_i)$ with errors $\Delta y_i$

 data = {{0.015, 34.1114},{0.0277, 35.705},{0.048948, 36.7316},{0.0651, 37.3067},{0.100915, 38.4567},{0.159, 39.4164},{0.248508, 40.2722},{0.455, 42.3239},{0.655, 42.3151},{0.75, 43.243},{0.84, 43.5143},{0.961, 44.2642},{1.188, 44.6076},{1.34, 45.0675},{1.414, 44.8038}}



that I want o fit with a complicated formula

$$Y(z) = 5\log \Bigl((1+z)\int_0^z \frac{dz'}{a+bz'+cz'^2}\Bigr)+25;$$

I used a NonLinearFit in the following way

model[a_?NumericQ, b_?NumericQ, c_?NumericQ, z_] := 5*Log[10, (1 + z)*
NIntegrate[1/(a + b*x + c*x^2), {x, 0, z}, PrecisionGoal -> 7,
AccuracyGoal -> 7]] + 25

fit1 =  NonlinearModelFit[data, model[a,b,c,z], {a, b, c}, z,


Mathematica was not able to make the fit, and the problem seems to be that the coefficients becomes complex. I decide then to use the analytical form of the Integral and then make the fit (to avoid the integration),

model2=5Log[10,(1+z)*(2ArcTan[(b+2*c*z)/(Sqrt[4*a*c-b^2])])/(Sqrt[\4*a*c-b^2])]+25;


but the same problem appears. There is any way to fix somehow the conditions over $\{a,b,c\}$ ??. The full set of data is about 580 points, its a Cosmological fot for Supernovae. I appreciate your help.

• I tired to use Assumptions over the $\{a,b,c\}$ to be Reals, but it! doesn't work – Alejandro Guarnizo Mar 25 '14 at 16:45
• Hi @b.gatessucks I tried to add more details – Alejandro Guarnizo Mar 25 '14 at 18:09
• You have (1/a + b x ..) should be 1/(a + b x .. ). also is the z^2 in the formula supposed to be (z_prime)^2 ? – george2079 Mar 27 '14 at 15:22
• Yes, you're right...z_prime should be also in the squared term! – Alejandro Guarnizo Mar 27 '14 at 15:24
• I see the parenthesis error got introduced in the edit history..best fix those things. – george2079 Mar 27 '14 at 15:26

Here is a solution -- I added a constraint and manually found some good initial values.

 g[z_] = Simplify[
5 Log[10,(1 + z ) Integrate[1/(a + b x + c x^2), {x, 0, z}] ] + 25 ,
Assumptions -> {z > 0, -b^2 + 4 a c > 0 }]


5 (5 + Log[10,-((2 (1 + z) (ArcTan[b/Sqrt[-b^2 + 4 a c]] - ArcTan[(b + 2 c z)/Sqrt[-b^2 + 4 a c]]))/Sqrt[-b^2 + 4 a c])])

 fit = NonlinearModelFit[data , {g[z],
-b^2 + 4 a c > 0  },  (*<-  here is your constraint that keeps things real *)
{{a, 0.000254}, {b, .000016}, {c, .000016}}, z]

Show[{ListPlot[data], Plot[ fit[z], {z, 0, Sqrt[2]}] }]

fit["BestFitParameters"]  (*without weight*)


{a -> 0.00022158364316055183, b -> 0.0001735648928280129, c -> 0.00003377302712833471}

The above works fine if I add Weights -> sigmadata as well (resulting in a very small change to the fit)

You may find this useful for "manually" finding some good initial values

   Manipulate[
Show[ {Plot[ g[z] /. {a -> am , b -> bm, c -> cm} , {z, 0, Sqrt[2]}],
ListPlot[data]}, PlotRange -> All] ,
{{am, .05}, 0, .1}, {{bm, .5}, 0, 1}, {{cm, .5}, 0, 1}]


## edit

looking closer it seems the fit is pushing up against the constraint:

b^2 - 4 a c /. fit["BestFitParameters"]


1.9057*10^-10

Going back to the formulation, consider the degenerate case:

 g[z_] = Simplify[
5 Log[10, (1 + z) Integrate[
1/(a + b x + (b^2/4/a ) x^2), {x, 0, z}]] + 25,
Assumptions -> {z > 0, a > 0, b > 0, c > 0}]


25 + (5 Log[(2 z (1 + z))/(2 a + b z)])/Log[10]

 NonlinearModelFit[data, {g[z], {b > 0, a > 0}}, {{a,0.000254}, {b,.000016}}, z]


also proves to give a nice looking fit to the data. You can also get yet another solution space by considering -b^2 + 4 a c > 0. You can combine the three cases using Piecewsise:

which can be then fed to NonlinearModelFit, with only the a>0,b>0,c>0 constraint.

• Thanks @Pillsy, but I couldn't find yet the solution. When I apply the fit for the real dataset, the fit (even is plotted fine) is complex!. I just evaluate the fit at some point and gives me values like '45 +0i'. There is other problem related with the Confidence Regions, and it seems is related to the initial values I put. Mathematica is not able to get the confidence regions if the problem is constrained. There is way to put the full data set here? – Alejandro Guarnizo Mar 28 '14 at 12:58
• @AlejandroGuarnizo george2709 answered your question; I just edited it so that NonlinearModelFit was spelled correctly (the l was capitalized in the original, which is perhaps the easiest typo in all of Mathematica). Maybe try using Chop` on the results of the fitting function? That discards tiny imaginary parts. – Pillsy Mar 28 '14 at 17:03