NonLinearFit with "complicated" integral model

Question

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


 sigmadata={0.215239,0.118114,0.175773,0.242628,0.121087,0.21481,0.152213,0.306006,0.188809,0.198402,0.471047,}

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, 
   Weights -> 1/(sigmadata)^2]; 

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.

Answer 1

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.