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

share|improve this question
    
you may check this question.mathematica.stackexchange.com/questions/5579/… –  s.s.o Jul 24 '13 at 19:23
    
Interesting, so Mathematica does not support output of chi-squared tests like this? It seems to me that most practicing scientists really care about chi-squared values... –  diracula Jul 24 '13 at 20:30
    
I agree that functionality should be there. But it's trivial enough to implement as well. –  s.s.o Jul 24 '13 at 21:11
    
Isn't NLM["EstimatedVariance"] what you are looking for? EstimatedVariance –  grbl Aug 3 '13 at 23:25
add comment

Know someone who can answer? Share a link to this question via email, Google+, Twitter, or Facebook.

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.