1
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ you may check this question.mathematica.stackexchange.com/questions/5579/… $\endgroup$ – s.s.o Jul 24 '13 at 19:23
  • $\begingroup$ 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... $\endgroup$ – diracula Jul 24 '13 at 20:30
  • $\begingroup$ I agree that functionality should be there. But it's trivial enough to implement as well. $\endgroup$ – s.s.o Jul 24 '13 at 21:11
  • $\begingroup$ Isn't NLM["EstimatedVariance"] what you are looking for? EstimatedVariance $\endgroup$ – grbl Aug 3 '13 at 23:25
3
$\begingroup$

This might be considered an extended comment but maybe it's an answer because I don't think you can get there from here: Your model is way too complex for the relationship between distance and the dependent variable and getting chi-squared statistics is the least of the concerns.

First, here's a log-log plot of one of the sets of data (and they all look pretty much the same):

Dat = Table[Data[x], {x, AU, 100 AU, AU}];
ListLogLogPlot[Dat]

Log log plot of data

The data looks pretty linear. A linear fit produces

lm = LinearModelFit[Log[Dat], x, x];
lm["BestFitParameters"]
(* {46.32985382377703262610476070259413913989`25.70992350462501, 
-1.99980476175735630477011151084940583553`25.70992350462501} *)

This is not surprising if we look at a series expansion of the log of the dependent variable given the log of the distance:

f = Log[GMsun/dist^2 + (GMjup dist)/(dist^2 + rJup^2)^(3/2) /. dist -> Exp[logDist]];
Series[f, {logDist, 0, 6}]

Series expansion The first two terms dominate and match very closely with the linear fit. In short, there's just not enough departure from a linear fit using the logs to estimate the parameters of the model of interest.

$\endgroup$
  • $\begingroup$ Thank you for responding. It has been a while since I asked this question, but if I remember correctly, the project I was working on was to determine that small deviation (alpha) from the linear fit. I was able to develop an independent chi-square test to solve this issue. $\endgroup$ – diracula Dec 10 '16 at 16:53
  • 1
    $\begingroup$ Glad you were successful. Note that hypothesis tests (chi-square tests and the like) have fallen out of favor over the last decade or so (usually for very good reasons) and there is more of an emphasis on estimation: standard errors, confidence intervals, etc. One of the reasons (that has always existed) is that the P-values produced by such tests can be very small due to large sample sizes as opposed to being due to the effect one is testing. Emphasizing estimation forces one to think about what size of effect is really important. Subject matter knowledge ought to be important. $\endgroup$ – JimB Dec 10 '16 at 17:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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