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$
4
  • $\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$
5
  • $\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
  • $\begingroup$ @JimB, you mention that hypothesis tests have fallen out of favour, would you be able to say what has replaced them if anything? I'm trying find the best way to characterise the goodness of a fit on a custom distribution and I was wondering what techniques are used today. Many thanks! $\endgroup$
    – alex
    Oct 19 '21 at 12:58
  • 1
    $\begingroup$ @alex Hypothesis tests are still "correct" but many times (in my opinion most of the time) doesn't give one an adequate decision making process in that an alternative hypothesis is rarely even considered. A significant result is the combination of the size of the departure from the null hypothesis AND the sample size. Without specifying the alternative hypothesis one can't separate out the effect of the sample size (which is irrelevant to the underlying issue). Getting an estimate of an effect and a measure of precision is much more useful. What should replace it? Thinking. $\endgroup$
    – JimB
    Oct 19 '21 at 15:01
  • $\begingroup$ ah I see. The way you phrased it sounded as if those tests have fallen out of favour and have already been replaced by a different (more favourable!) method. For someone who has a rather rudimentary understanding of statistics, it would be useful to keep up to date with the new techniques so if you think of something, please let me know :) $\endgroup$
    – alex
    Oct 19 '21 at 15:34

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.