I am using NonLinearModelFit for some curve fitting and I was wondering if NLM is able to output chi-squared/leastsquared statistics from the best-fit parameters and confidence intervals. From my understanding, NLM uses a least squares algorithm to find the best parameters, so shouldn't there be an associated chi-squared value with the fit?
To check that everything was working, I ran the test a thousand times and looked at the distribution of a certain parameter, alpha. Instead of finding a normal distribution however, I found a distribution with 2 peaks. To me, this indicates that: 1 the fitting function is messing up somehow, or 2, that it is choosing a local minimum and not searching for a better place. I have tried increasing the number of iterations so that the algorithm could possibly find a better local minimum, but that was not successful. So I currently think that the algorithm is getting caught up somewhere. So, I was hoping that I could check the chi-squared value for each iteration and see if some fits were better than others.
My code is below:
AU = 149597871000;
G = 6.67428*10^-11;
GMsun = 1.32712442099*10^20;
GMjup = GMsun/1047.348644;
dela = 10^-10;
rJup = 5.2 AU;
lambda = AU;
precision = 25;
alphas = {};
Data[dist_] := {SetPrecision[dist, precision],
SetPrecision[
GMsun/dist^2 + (GMjup dist)/(dist^2 + rJup^2)^(3/2) +
RandomReal[NormalDistribution[]] dela, precision]};
Model[dist_, alpha_, jupiter_, sun_, lambda_] :=
SetPrecision[(G sun)/dist^2 (1 + alpha Exp[-dist/lambda]) + (
G jupiter dist)/(dist^2 + rJup^2)^(3/2), precision];
Do[
Dat = Table[Data[x], {x, AU, 100 AU, AU}];
NLM = NonlinearModelFit[
Dat, {Model[dist, alpha, jupiter, sun, lam]}, {{alpha,
10^-7}, {jupiter, GMjup/G}, {sun, GMsun/G}, {lam, 20*AU}}, dist,
Tolerance -> 10^-50, AccuracyGoal -> precision,
WorkingPrecision -> precision, MaxIterations -> 1000];
realAlpha = NLM["ParameterTableEntries"][[1]][[1]];
realLambda = NLM["ParameterTableEntries"][[4]][[1]];
realJupiter = NLM["ParameterTableEntries"][[2]][[1]];
realSun = NLM["ParameterTableEntries"][[3]][[1]];
alphas = Append[alphas, Abs[realAlpha]];
, {i, 1000}]
Here the list alphas contains 1000 best-fit alphas from 1000 artificially created data sets (Note: this takes a while to run). The problem is that almost 20% of the time it outputs alpha ~ 10^-3, which is much too large to make sense.
Thank you!