0
$\begingroup$

I need some help to understand the proper meaning of the way Mathematica estimates interval of confidence. I have some data points, which are quite scattered around a sinusoidal behavior. For physical reasons, I expect the model to look like $$y(t)=p_0 + p_1\times \cos(\frac{2\pi}{365}t + \phi)$$

and I obtained with NonlinearModelFit $p_0 = 0.13$, $p_1=-0.07$ and $\phi = 0.04$.

Figure below show data points, the best fit (solid line), the SinglePredictionBands at 95% (blue band) and the MeanPredictionBands at 95% (yellow band).

enter image description here

My questions are the following

  1. What I would like to picture is : at day $D$, what is the interval in which I have 95% to find a value if I were to perform a measurement ? Does this corresponds to SinglePredictionBands ?
  2. If so, why isn't the interval closer to the fit value around, say, days 320-350 where data are much less spread then around day 175 ?
  3. If not, how can I get this interval ?
  4. What is MeanPredictionBands actually corresponding to ? Is it something like what the other reasonable fit would look like, considering errors on $p_0$, $p_1$ and $\phi$ ?

Here are the data and code (I took out some cosmetic considerations to keep the code a MWE)

data={7579/99136, 16595/297408, 3907/99136, 3463/74352, 19573/297408, 5173/99136, 9131/148704, 16591/297408, 2211/49568, 12529/297408, 7797/99136, 6413/74352, 17663/297408, 3065/49568, 5719/148704, 2179/99136, 3289/49568, 9209/148704, 2123/74352, 1561/37176, 13429/297408, 8147/148704, 8947/148704, 1283/24784, 2065/24784, 8729/297408, 3271/74352, 5131/99136, 1007/49568, 8063/148704, 2925/49568, 18787/297408, 10519/148704, 3043/74352, 4413/99136, 5893/99136, 3415/74352, 29549/297408, 13429/148704, 36449/297408, 8489/74352, 19159/148704, 29819/297408, 2079/24784, 14881/297408, 16255/297408, 8877/99136, 19727/297408, 9831/99136, 38677/297408, 4827/49568, 1531/24784, 18395/148704, 6995/99136, 8069/74352, 7773/99136, 6091/99136, 313/3098, 2591/24784, 14551/297408, 4129/49568, 15317/148704, 991/9294, 51949/297408, 4503/24784, 11695/74352, 24391/148704, 6227/49568, 37963/297408, 42785/297408, 15529/99136, 30245/297408, 28121/297408, 23369/297408, 15643/148704, 21503/148704, 38539/297408, 39257/297408, 10897/148704, 5197/99136, 2927/37176, 871/6196, 20681/297408, 7295/74352, 36553/297408, 14931/99136, 28621/297408, 483/6196, 28655/297408, 18679/148704, 12389/74352, 2555/18588, 36629/297408, 3661/37176, 61981/297408, 21657/99136, 33541/148704, 33883/148704, 2703/12392, 17883/99136, 13863/99136, 65093/297408, 32989/148704, 66809/297408, 312/1549, 32861/297408, 28741/297408, 8541/49568, 40487/297408, 18365/99136, 8797/37176, 55933/297408, 52207/297408, 19015/99136, 3039/24784, 5161/49568, 6923/74352, 18667/99136, 28489/148704, 48479/297408, 6269/99136, 12323/99136, 17809/148704, 17399/148704, 14327/99136, 29383/148704, 21167/99136, 1785/12392, 18259/99136, 2887/12392, 32461/148704, 20463/99136, 56687/297408, 21131/148704, 22415/148704, 1985/12392, 65333/297408, 22303/99136, 41113/297408, 14419/74352, 63817/297408, 8341/37176, 51179/297408, 54791/297408, 47617/297408, 4357/24784, 8975/37176, 23017/99136, 59795/297408, 56749/297408, 7027/37176, 50897/297408, 63629/297408, 689/3098, 12019/49568, 69047/297408, 67193/297408, 8365/37176, 33623/148704, 19763/99136, 41029/297408, 55657/297408, 1567/12392, 53557/297408, 44849/297408, 15471/99136, 44143/297408, 36695/148704, 56095/297408, 31663/148704, 69263/297408, 22227/99136, 15851/74352, 65879/297408, 76763/297408, 24707/99136, 71513/297408, 21187/99136, 37669/148704, 3033/12392, 24249/99136, 8591/37176, 29071/148704, 2053/9294, 22261/99136, 1295/6196, 8879/37176, 30175/148704, 1721/9294, 70601/297408, 24887/99136, 7895/37176, 64151/297408, 63383/297408, 64907/297408, 22439/99136, 64111/297408, 30889/148704, 18391/99136, 52169/297408, 56231/297408, 10065/49568, 15497/99136, 10793/49568, 53309/297408, 898/4647, 11/64, 16993/99136, 53837/297408, 15263/99136, 4309/24784, 16785/99136, 52249/297408, 36481/148704, 1393/6196, 23113/148704, 33845/148704, 64703/297408, 9643/49568, 30649/297408, 13103/99136, 9113/49568, 308/1549, 20623/99136, 2287/18588, 34633/297408, 5731/37176, 19169/99136, 9115/49568, 54299/297408, 10417/49568, 10379/49568, 5643/24784, 7031/37176, 16321/148704, 19519/148704, 19803/99136, 10235/49568, 15435/99136, 18067/99136, 15535/74352, 20457/99136, 10789/74352, 14303/99136, 52535/297408, 43787/297408, 15151/99136, 7941/49568, 57721/297408, 64009/297408, 19861/99136, 20679/99136, 52153/297408, 8513/49568, 26347/297408, 13007/148704, 649/4647, 28471/297408, 30317/297408, 2111/18588, 11897/74352, 14599/99136, 19095/99136, 18805/99136, 30769/297408, 46211/297408, 2063/12392, 47389/297408, 53105/297408, 53989/297408, 9119/49568, 49159/297408, 25487/148704, 8903/49568, 12585/99136, 13537/148704, 3/32, 11993/148704, 28991/297408, 6097/49568, 38365/297408, 6047/37176, 3823/24784, 14129/99136, 5693/74352, 8215/99136, 12941/99136, 15551/148704, 12765/99136, 1109/9294, 8543/148704, 743/6196, 33265/297408, 15143/148704, 3035/49568, 32135/297408, 23855/297408, 29755/297408, 32573/297408, 16409/297408, 4943/74352, 2671/24784, 10009/74352, 43249/297408, 449/3098, 31453/297408, 5473/99136, 2029/24784, 4935/49568, 2457/24784, 18385/148704, 34853/297408, 445/4647, 29119/297408, 14183/148704, 3239/37176, 5523/99136, 28733/297408, 2695/24784, 21623/297408, 19979/297408, 24235/297408, 17333/297408, 11483/297408, 4775/74352, 27811/297408, 24443/297408, 1665/24784, 889/24784, 17441/297408, 5797/74352, 4295/49568, 2703/49568, 11543/148704, 22535/297408, 8293/99136, 10471/148704, 8299/148704, 9157/148704, 1571/24784, 1523/18588, 635/12392, 8731/99136, 2105/24784, 9941/148704, 6169/74352, 19607/297408, 19691/297408, 973/24784, 7411/148704, 8579/148704, 5125/99136, 20125/297408, 767/18588, 6197/99136, 10211/148704, 15901/297408, 5449/99136, 7589/99136, 4193/49568, 21011/297408, 8257/148704, 5231/99136, 12095/148704, 5855/148704};

model = P0 + P1*Cos[(2 \[Pi])/365*t + \[Phi]];
fit = NonlinearModelFit[data, 
model, {{P0, 1}, {P1, 500}, {\[Phi], \[Pi]/2}}, t]
band95S[t_] = fit["SinglePredictionBands", ConfidenceLevel -> .95];
band95M[t_] = fit["MeanPredictionBands", ConfidenceLevel -> .95];
Show[
    Plot[{{fitSol[x], band95M[x], band95S[x]}}, {x, 0, 365},Filling -> {{2 -> {1}}, {3 -> {1}}}], 
    ListPlot[data, PlotStyle ->Directive[PointSize[Medium]]], 
    PlotRange -> {{0, 365}, {0, 0.28}}, Frame -> True]
Improve this question
$\endgroup$
2
  • $\begingroup$ Please provide the data and code you have so far. $\endgroup$
    – Feyre
    Feb 12, 2017 at 14:41
  • $\begingroup$ @Feyre : done. Hope this helps. $\endgroup$
    – DSuchet
    Feb 13, 2017 at 1:21

1 Answer 1

Reset to default
5
$\begingroup$

I would argue that all of your questions apply to any statistical regression package and the answers are all about statistical regression analysis concepts rather than Mathematica.

With NonlinearModelFit you are fitting the model

$$y(t)=p_0+p_1 \cos{(2 \pi t/365+\phi)}+\epsilon_t$$

with independent and identically distributed errors where $\epsilon_t \sim N(0,\sigma^2)$ for $t=1,2,\ldots,365$.

  1. What I would like to picture is : at day DD, what is the interval in which I have 95% to find a value if I were to perform a measurement ? Does this corresponds to SinglePredictionBands ?

95% of the time future observations will be between

$$p_0+p_1 \cos{(2\pi t/365+\phi)} \pm 1.96\sigma$$

given the value of $t$ assuming the above model is the actual model generating the data. This is a non-unique interval where 95% of the future observations will occur. However, we don't know the values of the parameters and need to estimate those parameters. The SinglePredictionBands account for the fact that all of the paramters (incliding the error variance) are estimated. If one would repeatedly sample the same 365 days, then 95% of these SinglePredictionBands at a value of $t$ would contain a future observation.

  1. If so, why isn't the interval closer to the fit value around, say, days 320-350 where data are much less spread then around day 175 ?

You and NonlinearModelFit have assumed that the variance is constant for all values of $t$. If the data does not appear to support that assumption, you need to fit a model that accounts for the varying of the error variance. (Error variances are parameters, too!)

  1. If not, how can I get this interval ?

You would need to model the error variance across days. NonlinearModelFit can do some kinds of such models using the Weights option or through a transformation of the response. But first you need to determine if you really do need to model the variance. Residual plots (such as the square root of the absolute residual values vs. the predicted value or vs. the day) might be informative.

  1. What is MeanPredictionBands actually corresponding to ? Is it something like what the other reasonable fit would look like, considering errors on $p_0$, $p_1$, and $\phi$ ?

What Mathematica calls MeanPredictionBands are 95% confidence intervals for

$$p_0+p_1 \cos{(2\pi t/365+\phi)}$$

Improve this answer
$\endgroup$

Your Answer

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

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