2
$\begingroup$

In the realm of TimeSeries of data with Circular Statistics, NonlinearModelFit (NLM) correctly generates reasonable Standard Error & confidence intervals for phase parameters (phaseCIs) with x-values that are integer offsets from a reference time (TRef). However, when AbsoluteTime values are utilized, resulting NLM Standard Error for phaseCIs are nonsensical. NOTE: I'm aware of FittedModel::constr warning, but still...

Questions:

  1. Why does NLM ignore parameter constraints Abs[phi] <= 2Pi when AbsoluteTime units are utilized?
  2. Is there an additional NLM tuning parameter to ensure Phi never steps outside its constrained param value?
  3. Is there a two-pass strategy to get reasonable phaseCIs (e.g., estimate BestFitParameters in one pass with a specific NLM Method, and then run a second NLM passing with a Method that explicitly applies parameter constraints)?
  4. (More generally) What other issues does NLM have with respect to TimeSeries data with Circular Statistical models?
(* Globals *)
{dayAsSeconds, lunarDayAsSeconds} = 
  QuantityMagnitude @  
     UnitConvert[Quantity[#, "Hours"], "Seconds"] & /@ {24, 24.8};
TRef = DateObject@ {2021, 1, 12, 0, 0} ;
deltaT = 
  QuantityMagnitude @  UnitConvert[Quantity[30, "Minutes"], "Seconds"];
(* Define xValues, one with TRef offset, the other AbsoluteTime *)
xAsOffsets =  Range[5 dayAsSeconds, 20 dayAsSeconds , deltaT ];
xAsAbsoluteTime = 
  AbsoluteTime[TRef + Quantity[#, "Seconds"] ] & /@ xAsOffsets;

NOTE: At this point, the only difference between xAsOffsets and xAsAbsoluteTime is their IntegerLength[] size: the former are offset values from a TRef {0, 3600, ...} with IntegerLength[] <= 7 whereas the latter are AbsoluteTime values {3819830400, 3819832200, ...} whose IntegerLength[] == 10. Both are monotonically increasing by 1800.

(* Model *)
model = M + A1 Cos[2 Pi (t / tau1) + phi1] +   
   A2 Cos[2 Pi (t / tau2) + phi2] ;
modelSeeds =  { M -> 100, A1 -> 10, tau1 -> dayAsSeconds, 
   phi1 -> - Pi/4, A2 -> 2, tau2 -> lunarDayAsSeconds, 
   phi2 -> - Pi/4};
SeedRandom[1];
rnd[dummy_] := 5 RandomReal [{-.5, .5}];
(* Generate xyPairs *)
xyAsOffset = 
  Table[{t, (model /. modelSeeds) + rnd[t]} , {t, xAsOffsets }];
xyAsAbsoluteTime = 
  Table[{t, (model /. modelSeeds) + rnd[t]} , {t, xAsAbsoluteTime }];
(* Prepare NLM wrapper *)
nlmParms = { {M}, {A1}, {tau1, dayAsSeconds}, {phi1, Pi}, {A2}, {tau2,
     lunarDayAsSeconds}, {phi2, Pi}};
nlmConstraint = 
  Abs[phi1] <= 2 Pi && Abs[phi2] <= 2 Pi && A1 > 0 && A2 > 0;
nlm[data_, opts : OptionsPattern[NonlinearModelFit]] := 
 nlm[data, opts] = 
  NonlinearModelFit[data, {model, nlmConstraint}, nlmParms, t, opts];

Visually verify both computed dataset Y-values are in the same phase

ListPlot[{xyAsOffset[[All, 2]], xyAsAbsoluteTime[[All, 2]]}, 
 Joined -> True]

Confirming both datasets are in phase

Visually compare phi1 and phi2 Standard Error and Confidence Intervals to see the crazy phaseCIs for xyAsAbsoluteTime.

compareCIResults[opts : OptionsPattern[NonlinearModelFit]] := 
 TableForm[{#["ParameterConfidenceIntervalTable"] & /@ {nlm[
      xyAsOffset, opts], nlm[xyAsAbsoluteTime, opts]}}, 
  TableHeadings -> {{opts}, {"xAsOffsetsFromTRef", 
     "xAsAbsoluteTime"}}];
compareCIResults[]

Comparison of ParameterConfidenceIntervalTable

$\endgroup$
1
  • $\begingroup$ You might get more help if you removed some of the unnecessary information. For example replacing xAsOffset and xAsAbsoluteTime as x1 and x2 with x1 = 432000 + 1800 Range[0, 720]; and x2 =x1 + 3819398400;. The issue you raise is legit. It's just clouded by a lot of unnecessary information. And I don't think this qualifies as "circular statistics" or whatever a "circular time series" is just because cosines are involved. It's simply a regression with an additive error. $\endgroup$
    – JimB
    Dec 18, 2021 at 5:29

1 Answer 1

3
$\begingroup$

This (for me) has been challenging but in essence, given the data and model, the parameters phi1 and tau1 are almost perfectly correlated (as are phi2 and tau2) which means one can set phi1 and phi2 to any value and you'll get an almost identical fit. So the "crazy" standard errors are due to the overparameterization of the model rather than any error in NonlinearModelFit.

Consider an example that exactly matches your data generation and model.

(* Specify model and parameter values of interest *)
model = m + a1 Cos[2 Pi t/(10^4 tau1) + phi1 ] + 
   a2 Cos[2 Pi t/(10^4 tau2) + phi2 ];
parameters = {m -> 100, a1 -> 10, tau1 -> 8.64, a2 -> 2, 
   tau2 -> 8.928, phi1 -> -π/4, phi2 -> -π/4};

Note that I've scaled tau1 and tau2 so that all parameters are nearly the same order of magnitude which is generally good practice and tends to give one more numerically stable results.

(* Specify predictor values *)
n = 701;
x = 381830400 + 1800 Range[0, n];

(* Generate response *)
SeedRandom[12345];
error = RandomVariate[UniformDistribution[{-2.5, 2.5}], n + 1];
y = (model /. parameters /. t -> x) + error;
xy = Transpose[{x, y}];

(* Fit regression *)
restrictions = {a1 > 0, a2 > 0, tau1 > 0, tau2 > 0};
(* Model fitting all parameters *)
nlm1 = NonlinearModelFit[xy, {model, restrictions}, 
  {{m, 100}, {a1, 10}, {tau1, 8.64}, {phi1, -Pi/4},
     {a2, 2}, {tau2, 8.928}, {phi2, -Pi/4}}, t] // Quiet;
(* Model fixing phi1 and phi2 - doesn't matter what values are used *)   
nlm2 = NonlinearModelFit[xy, {model /. {phi1 -> 1, phi2 -> 1}, restrictions},
    {{m, 100}, {a1, 10}, {tau1, 8.64}, {a2, 2}, {tau2, 8.928}}, t] // Quiet;

The fits are nearly identical.

(* Show results *)
Show[Plot[{nlm1[t], nlm2[t]}, {t, Min[x], Max[x]}, 
  PlotStyle -> {{LightGray, Thickness[0.01]}, {Red, Thickness[0.0005]}}],
  ListPlot[xy], AspectRatio -> 1/4, ImageSize -> Full]

Data and fits

(One fit is in light gray with a wider line and the other fit is in red.)

To justify the claim that phi1 and tau1 (and phi2 and tau2) are nearly perfectly correlated consider the correlation matrix from nlm1 above (where all parameters were fit):

NumberForm[nlm1["CorrelationMatrix"] // MatrixForm, 3]

Correlation matrix

There are other highly correlated estimators. In fact, all parameters except m are highly correlated with each other which can indicate numerical instability in the estimates of those parameters.

$\endgroup$
2
  • $\begingroup$ JimB -- Thank you for simplifying this model -- and educating me on the correct strategy for generating the random error correctly! That said, off-line I'll re-work the example and if I may, send directly to you so that I'm clear about your point with the over-parameterization. $\endgroup$ Dec 31, 2021 at 3:05
  • $\begingroup$ JimB -- ALSO thanks for pointing out the tip of rescaling Tau1 and Tau2 -- more for me to fold into on-going tests here as I explore your response further here. $\endgroup$ Jan 1, 2022 at 21:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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