0
$\begingroup$

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
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Bob Hanlon
    May 1 at 15:11
  • $\begingroup$ 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... $\endgroup$ May 1 at 15:21
  • $\begingroup$ Rather than drop the WorkingPrecision, either use exact values or arbitrary-precision numbers throughout. $\endgroup$
    – Bob Hanlon
    May 1 at 15:49

1 Answer 1

1
$\begingroup$

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

Parameter correlation matrix

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]]]}]]

Plot of data and predicted response

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

Warning message about invalid non-numeric values

$\endgroup$
5
  • $\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$ May 1 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
    May 1 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$ May 1 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
    May 1 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$ May 2 at 0:37

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.