Interpreting Confidence Intervals From NonLinearModelFit

Question

I recently applied a model to some example data using nonlinearmodelfit and, while the fit was seemingly good (R^2=0.987) and the returned parameters were mostly in line with expectations without constraining the model, the confidence intervals were all over the place. Some returned parameters had reasonable confidence intervals (m1,m2,ecc,t01) but the rest had these ridiculous error bars that seem like trouble. Is this evidence that there is something wrong with the model/parameters or is there something with nonlinearmodelfit that needs adjusting? The code is down below. Apologies that there is a lot going on in it but it should only take a few seconds to run.

ClearAll["Global`*"]
lc = {{0.`, 0.05945360824742174`}, {0.599398`, 
    0.06045360824742119`}, {1.198796`, 
    0.042453608247420505`}, {1.7981939999999998`, 
    0.038453608247420945`}, {2.397592`, 
    0.021453608247421485`}, {2.9969900000000003`, 
    0.030453608247421826`}, {3.5963879999999997`, 
    0.007453608247420362`}, {4.195786`, 
    0.0244536082474216`}, {4.795184`, 
    0.004453608247420249`}, {5.394582`, -0.011546391752579765`}, \
{5.9939800000000005`, -0.017546391752578216`}, {6.5933779999999995`, \
-0.04054639175257968`}, {7.192775999999999`, -0.04054639175257968`}, \
{7.792174`, -0.060546391752579254`}, {8.391572`, \
-0.04754639175257935`}, {8.990969999999999`, -0.06354639175257937`}, \
{9.590368`, -0.06754639175257893`}, {10.189766`, \
-0.07754639175257871`}, {10.789164`, -0.09354639175257873`}, \
{11.388562`, -0.0825463917525795`}, {11.987960000000001`, \
-0.08054639175257883`}, {12.587357999999998`, -0.08154639175257827`}, \
{13.186755999999999`, -0.09354639175257873`}, {13.786154`, \
-0.07354639175257915`}, {14.385551999999999`, -0.08654639175257905`}, \
{14.98495`, -0.07854639175257816`}, {15.584348`, \
-0.06254639175257815`}, {16.183746`, -0.06754639175257893`}, \
{16.783144`, -0.0615463917525787`}, {17.382541999999997`, \
-0.05054639175257947`}, {17.981939999999998`, -0.04654639175257813`}, \
{18.581338`, -0.024546391752579666`}, {19.180736`, \
-0.014546391752578103`}, {19.780134`, -0.010546391752578543`}, \
{20.379532`, 0.005453608247421471`}, {20.97893`, 
    0.015453608247421258`}, {21.578328`, 
    0.028453608247421158`}, {22.177726`, 
    0.04345360824742173`}, {22.777124`, 
    0.04845360824742073`}, {23.376522`, 
    0.06145360824742063`}, {23.975920000000002`, 
    0.07245360824742164`}, {24.575318`, 
    0.09345360824742066`}, {25.174715999999997`, 
    0.10045360824742033`}, {25.774113999999997`, 
    0.10945360824742068`}, {26.373511999999998`, 
    0.11245360824742079`}, {26.97291`, 
    0.10945360824742068`}, {27.572308`, 
    0.1104536082474219`}, {28.171705999999997`, 
    0.11145360824742134`}, {28.771103999999998`, 
    0.11145360824742134`}, {29.370502`, 
    0.10545360824742112`}, {30.569298`, 
    0.08245360824742143`}, {31.168696`, 
    0.0894536082474211`}, {31.768094`, 
    0.05945360824742174`}, {32.367492`, 
    0.06645360824742141`}, {32.96689`, 
    0.04545360824742062`}, {33.566288`, 
    0.03245360824742072`}, {34.165685999999994`, 
    0.019453608247420817`}, {34.765083999999995`, 
    0.003453608247420803`}, {35.364481999999995`, \
-0.010546391752578543`}, {35.963879999999996`, \
-0.023546391752578444`}, {36.563278`, -0.03154639175257934`}, \
{37.162676`, -0.032546391752578785`}, {37.762074`, \
-0.04754639175257935`}, {38.361472`, -0.0615463917525787`}, \
{39.560268`, -0.07354639175257915`}, {40.159666`, \
-0.07954639175257938`}, {40.759064`, -0.09154639175257984`}, \
{41.358461999999996`, -0.0825463917525795`}, {41.95786`, \
-0.10454639175257974`}, {43.156656`, -0.10154639175257962`}, \
{43.756054`, -0.0745463917525786`}, {44.355452`, \
-0.07754639175257871`}, {44.95485`, -0.08554639175257961`}, \
{45.554248`, -0.07754639175257871`}, {46.153646`, \
-0.07254639175257971`}, {46.753044`, -0.045546391752578685`}, \
{47.352442`, -0.05154639175257891`}, {47.951840000000004`, \
-0.03354639175257823`}, {48.551238000000005`, \
-0.028546391752579225`}, {49.150636`, -0.014546391752578103`}, \
{49.750034`, -0.011546391752579765`}, {50.34943199999999`, 
    0.004453608247420249`}, {50.948829999999994`, 
    0.008453608247421585`}, {51.548227999999995`, 
    0.019453608247420817`}, {52.147625999999995`, 
    0.021453608247421485`}, {52.747023999999996`, 
    0.03545360824742083`}, {53.346422`, 
    0.04645360824742184`}, {53.94582`, 
    0.05745360824742107`}, {54.545218`, 
    0.06145360824742063`}, {55.144616`, 
    0.06445360824742075`}, {55.744014`, 
    0.07045360824742097`}, {56.343411999999994`, 
    0.06745360824742086`}, {56.942809999999994`, 
    0.07345360824742109`}, {57.542207999999995`, 
    0.0764536082474212`}, {58.141605999999996`, 
    0.07745360824742065`}, {58.741004`, 
    0.0684536082474203`}, {59.340402`, 0.06945360824742153`}};
error = Table[0.006729440487569607, {x, 0, 96}];
G = 6.6743*10^-11;
c = 2.99792458*10^8;
Periods = 59.9394;
P = Periods*86400;
a = Power[(P^2*G*(m1 + m2))/(4*\[Pi]^2), (3)^-1]
\[Mu] = (m1 m2)/(m1 + m2)
aph = a (1 + ecc);
x0 = (m2*aph)/(m1 + m2)
y0 = 0;
X0 = -((m1*aph)/(m1 + m2))
Y0 = 0;
Vel1 = ((G \[Mu] )/a  m2/m1 (1 - ecc)/(1 + ecc) )^(
 1/2)(*initial velocity star 1*)
Vel2 = -((G \[Mu] )/a  m2/m1 (1 - ecc)/(1 + ecc) )^(1/2) m1/
  m2(*initial velocity star 2*)


system1 = {
   x1''[t] == -((G m1 m2  (x1[t] - x2[t]))/(
     m1 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))),
   y1''[t] == -((G m1 m2  (y1[t] - y2[t]))/(
     m1 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))),
   x2''[t] == (G m1 m2  (x1[t] - x2[t]))/(
    m2 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2)),
   y2''[t] == (G m1 m2  (y1[t] - y2[t]))/(
    m2 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))};

initials = {
   x1[0] == x0,
   x1'[0] == 0,
   y1[0] == 0,
   y1'[0] == Vel1,
   x2[0] == X0,
   x2'[0] == 0,
   y2[0] == 0,
   y2'[0] == Vel2};

{xbh1, ybh1, xbh2, ybh2} = 
 ParametricNDSolve[{system1, initials}, {x1, y1, x2, y2}, {t, 0, 
   P}, {m1, m2, ecc}]
bh1 = {(xbh1[t*86400]/(1.5*10^11)), (ybh1[t*86400]/(1.5*10^11))};
bh2 = {(xbh2[t*86400]/(1.5*10^11)), (ybh2[t*86400]/(1.5*10^11))};
G = 6.6743*10^-11; (*Gravitational Constant*)
M = m1; (*Mass of BH ()*)
m = m2; (*Mass of Source BH (kg)*)
c = 2.99792458*10^8; (*speed of light*)
xp1 = x1[m1, m2, ecc][(t + t01)*86400] /. xbh1;
yp1 = y1[m1, m2, ecc][(t + t01)*86400] /. ybh1;
xp2 = x2[m1, m2, ecc][(t + t01)*86400] /. xbh2;
yp2 = y2[m1, m2, ecc][(t + t01)*86400] /. ybh2;
vxp2 = D[xp1, t]/86400;
vyp2 = D[xp2, t]/86400;
\[Omega] = 1.5428711;
\[Alpha] = 1;
g = 0.46;
\[CapitalGamma] = 0.64;
sev = (0.15 (15 + \[CapitalGamma]) (1 + g))/(3 - \[CapitalGamma]) m1/
   m2 (rstar (6.957*10^8)/\[Sqrt]((xp2 - xp1)^2 + (yp2 - 
           yp1)^2))^3 Cos[
    2 ((2 \[Pi])/P*((t + t01)*86400) + \[Pi] + \[Omega])] Sin[\[Pi]/
     2 - \[Phi]]^2;

Delf = ((3 - (-\[Alpha])) (((1/
         c (-(vxp2*Sin[\[Omega]] Sin[\[Pi]/2 - \[Phi]] + 
             vyp2*Cos[\[Omega]] Sin[\[Pi]/2 - \[Phi]]))))*
      Sin[\[Pi]/2 - \[Phi]])) ;
myfit = NonlinearModelFit[lc, 
  Re[0 - 2.5 Log10[((1 + (Delf + sev)))]], {{m1, 
    2.9 (1.988*10^30)}, {m2, 0.74 (1.988*10^30)}, {rstar, 25}, {ecc, 
    0.02}, {t01, 2.9}, {\[Phi], 0.05}}, {t}, Weights -> 1/error^2, 
  ConfidenceLevel -> 0.68, 
  VarianceEstimatorFunction -> (Mean[Abs[#]] &), Tolerance -> 10^-50]
myfit["ParameterTable"]
myfit["ParameterConfidenceIntervalTable"]
Grid[Transpose[{#, myfit[#]} &[{"AdjustedRSquared", "AIC", "BIC", 
    "RSquared"}]], Alignment -> Left]
myfit["CorrelationMatrix"] // Quiet // MatrixForm

Please check the code that you posted. It produces numerous error messages and hangs longer than I'm willing to wait. Also, it makes no sense to ask for WorkingPrecision -> 17 in the NonlinearModelFit with machine precision input. — Bob Hanlon, Commented May 1, 2022 at 15:11
Hmmm interesting. Fair enough with the working precision (I just edited it out in the code above), though now the errorbars are large for everything (an interesting development). However, I copied both the original code and this new one from here and pasted it directly into mathematica and found it only took about 7 seconds to run. Not sure why it's hanging for you... — sorabella91, Commented May 1, 2022 at 15:21
Rather than drop the WorkingPrecision, either use exact values or arbitrary-precision numbers throughout. — Bob Hanlon, Commented May 1, 2022 at 15:49

JimB · Accepted Answer · 2022-05-01 16:31:32Z

1

Your model is likely over-parameterized in that estimators for $m_1$ and $m_2$ as are $\phi$ and $\alpha$ perfectly negatively correlated. This also typically results in large confidence intervals for those parameters.

myfit["CorrelationMatrix"] // Quiet // MatrixForm

This occurs if the data just doesn't support separating those pairs of parameters or if those parameters are under sum restriction such as $m_1+m_2=1$ or if $m_1 + m_2$ always occur together and are never found separately in the model definition.

From the "ParameterTable" it appears that only 2 (ecc and t01) of the 9 parameters are statistically significant which is another strong suggestion that the model is over-parameterized.

If interpreting the coefficients is your objective, then you've got trouble.

However, if prediction of the response variable is your objection, then an overparameterized model isn't so much trouble. A plot of the data and predicted response is

Show[ListPlot[lc], Plot[myfit[z], {z, Min[lc[[All, 1]]], Max[lc[[All, 1]]]}]]

There are still lots of warning messages that should be looked into such as:

answered May 1, 2022 at 16:31

JimB

42.9k3 gold badges51 silver badges108 bronze badges

$\begingroup$ Thank you for the in-depth reply! Yes I was worried about over-parameterization so thank you for confirming that suspicion. That warning message stems from trying to get the parameters on or near the same order of magnitude (I was told vastly different order of magnitudes could lead to strange results). Haven't figured out yet why that message appears and why it doesn't seem to stop the code from running. $\endgroup$
– sorabella91
Commented May 1, 2022 at 16:47
$\begingroup$ Getting the parameter estimates on the same order of magnitude is a good thing. But that warning message says that there is a non-numeric value. I would take that warning message at its word and not think it was about scaling coefficients. $\endgroup$
– JimB
Commented May 1, 2022 at 16:51
$\begingroup$ So, taking your suggestions, I edited the code above. I got rid of the order of magnitude stuff and the majority of the error messages went with it. I decided to treat a few of the variables as knowns instead of having the model fit to them (there are now only 6 parameters to be estimated instead of 9). Of the 6, five return estimates with significant p values (a better sign than before). However, with m1 and m2 on a significantly different order of magnitude, their standard error is practically zero (10^12 is not particularly comparable to 10^30). $\endgroup$
– sorabella91
Commented May 1, 2022 at 18:12
$\begingroup$ Hopefully someone else can help here. I just don't know the physics to be able to debug the issue(s). $\endgroup$
– JimB
Commented May 1, 2022 at 21:45
$\begingroup$ No problem. Thank you for the help though! It is strange as this model has worked for another data set without much issue. Hopefully I'll get to the bottom of it soon. $\endgroup$
– sorabella91
Commented May 2, 2022 at 0:37

Add a comment |

Stack Exchange Network

Interpreting Confidence Intervals From NonLinearModelFit

1 Answer 1

Your Answer

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

Hot Network Questions

Interpreting Confidence Intervals From NonLinearModelFit

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

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

Related

Hot Network Questions

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