Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
I tired to use Assumptions over the $\{a,b,c\}$ to be Reals, but it! doesn't work –  Alejandro Guarnizo Mar 25 at 16:45
    
Hi @b.gatessucks I tried to add more details –  Alejandro Guarnizo Mar 25 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 at 15:22
    
Yes, you're right...z_prime should be also in the squared term! –  Alejandro Guarnizo Mar 27 at 15:24
    
I see the parenthesis error got introduced in the edit history..best fix those things. –  george2079 Mar 27 at 15:26
show 6 more comments

1 Answer 1

up vote 4 down vote accepted

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

enter image description here

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:

enter image description here

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

share|improve this answer
    
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 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 at 17:03
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.