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.

First of all, a bit of background, I am using Mathematica to fit angular distributions of gamma rays, so I apologize for some of the serious messiness in the formulas.

The data is simple I have three lists. The first is dataWithErrors = {{x1,y1,dy1}, {x2, y2, dy2}, ... }. The second, rawData = {{x1,y1}, {x2,y2}, ... } is what actually goes to the model fit function. The third is ptWeights={1/dy1^2, 1/dy2^2, ... } Now the code:

P[m_, J1_, sig_]:=E^(-m^2/(2*sig^2))/Sum[E^(-mp^2/(2*sig^2)),{mp,-J1,J1}]
B[k_,J1_,sig_]:=((2*J1+1)^(1/2))*Sum[((-1)^(J1+m))*ClebschGordan[{J1,-m},{J1,m},{k,0}]*P[m,J1,sig],{m,-J1,J1}]
F[k_,L1_,L2_,J1_,J2_]:=((-1)^(J1+J2-1))*(((2*L1+2)*(2*L2+1)*(2*J2+1)*(2*k+1))^(1/2))*ThreeJSymbol[{L1,1},{L2,-1},{k,0}]*SixJSymbol[{L1,L2,k},{J2,J2,J1}]
A[k_,delt_,L1_,L2_,J1_,J2_]:=(F[k,L1,L2,J2,J1]+2*delt*F[k,L1,L2,J2,J1]+(delt^2)*F[k,L1,L2,J2,J1])/(1+delt^2)
ForTransition[k_,delt_,sig_]:=A[k,delt,1,2,17/2,15/2]*B[k,17/2,sig]
(*Constants to account for the detectors not being points*)
saCorr2 = 0.987751
saCorr4 = 0.959561

Now to perform the fit I use the following command:

TransitionFit = NonlinearModelFit[rawData,a0(1 + saCorr2*ForTransition[ 2, delt, sig]*LegendreP[2, Cos[th] ] + saCorr4*ForTransition[ 4, delt, sig]*LegendreP[4, Cos[th] ],{{a0,12544.9},{delt,-1.2395},{sig,1.983}}, th,Weights -> ptWeights, VarianceEstimatorFunction -> (1 &)]

I calculate the reduced chi^2 using:

(Plus@@(((#[[2]]-TransitionFit[#[[1]]])^2/(#[[3]]^2))&/@dataWithErrors)/(16-3)

Dividing by 16-3 because I have 16 points and 3 parameters.

Here is the problem, if I run the fit with sig as a free parameter I get a0=12567.8, delt=-1.28286, sig=1.94335, and a reduced chi^2 of 2.912.

All of which looks somewhat reasonable. But if I fix sig = 1.983 (and calculate reduced chi^2 exactly the same way) I get a0=12576.4, delt=-1.25549, and a reduced ch^2 of 2.901.

If I also fix delt = -1.23947... I get a0=12568.1 and a reduced chi^2 of 2.89344.

Finally if I simply calculate the chi^2 with a0 = 12544.9 and delt and sig as above I get a reduced chi^2 of 2.88363.

As you can see in the NonlinearModelFit call I give it the ideal values, and it still is not giving those values

I realize that these are small differences in the chi^2, but I don't understand why mathematica is doing this. Why is this happening? Are there any suggestions to fix it? Things I can do to improve how my code?

share|improve this question
add comment

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.