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

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

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