# NonlinearModelFit is not giving the fit with minimum chi^2

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?

-