Unexpected behavior of confidence bands for data presenting two regions with uneven noise

Question

Let's generate some noisy data with uneven noise.

    data = Table[{x,Exp[-(x-2)^2] + Exp[-(x+2)^2]*RandomReal[{0.5, 1.5}]},
           {x, RandomReal[{-4, 4}, 500]}];

Let's fit it with a nonlinear model

    fit = NonlinearModelFit[data, b Exp[-(x-a)^2] + c Exp[-(x+a)^2],
          {a, b, c}, x]

from which we extract the σ and 2σ prediction bands for single observations:

    {bands1[x_], bands2[x_]}=Table[fit["SinglePredictionBands",ConfidenceLevel -> cl],
      {cl, {.683, .954}}];

The bands seem not to take into account the unevenness of the noise:

    Show[ListPlot[data, PlotStyle -> {Opacity[0.5]}], Plot[{fit[x], bands1[x], bands2[x]},
         {x, -4, 4},Filling -> {2 -> {{1}, Directive[{Green, Opacity[0.25]}]}, 3 -> {{2}, Directive[{Yellow, Opacity[0.25]}]}},
         PlotStyle -> {Directive[Thickness[0.0075], Red],Black, Black}]]

Am I misinterpreting confidence bands, or am I calculating them in the wrong way?

Answer 1

I believe there are some big misunderstandings in the question...

First: usually error terms are additive, not multiplicative.

There are the so-called multiplicative error models (MEM), but I don't think it's the case here, because usually in MEM the multiplicative error is time-dependent (i.e., $\epsilon_{t} \in [0, \infty)$), while in the OP's question the error term is withdrawn from a fixed interval (i.e, $\epsilon \in [0.5, 1.5]$). A good introduction to MEM can be found here.

Second: This said, it seems that we cannot find confidence intervals for this "model" because it is actually possible to find certainty intervals (sorry if this term doesn't exist, maybe I've just created it...).

Consider, for instance, the original model:

$$ f(x)=\dfrac{1}{e^{(x-2)^2}} + \dfrac{1}{e^{(x+2)^2}}\epsilon, $$ s.t. $\epsilon \in [0.5, 1.5]$.

If we make a small change to the original model it will be easier to understand how to draw the "certainty intervals". Consider the following modification:

$$ f(x)=\dfrac{1}{e^{(x-2)^2}}\epsilon_{1} + \dfrac{1}{e^{(x+2)^2}}\epsilon_{2}. $$

If we set $\epsilon_{1} \in [1, 1]$ and $\epsilon_{2} \in [0.5, 1.5]$ we have the original model.

Now we can show that the certainty intervals are dependent upon the error terms ($\epsilon_{1}$ and $\epsilon_{2}$) intervals.

Manipulate[
  Module[{f, g, h, y}, 
  {f[y_] := error1*Exp[-(y - 2)^2], 
  g[y_] := (error1 + error2)/2*Exp[-(y - 2)^2], 
  h[y_] := error2*Exp[-(y - 2)^2], 
  G1 = Plot[{f[y], g[y], h[y]}, {y, -4, 4}, Filling -> {1 -> {3}}, 
    PlotRange -> {{-4, 4}, {0, 2}}]}];
  Module[{f, g, h, y}, {f[y_] := error3*Exp[-(y + 2)^2], 
  g[y_] := (error3 + error4)/2*Exp[-(y + 2)^2], 
  h[y_] := error4*Exp[-(y + 2)^2], 
  G2 = Plot[{f[y], g[y], h[y]}, {y, -4, 4}, Filling -> {1 -> {3}}, 
    PlotRange -> {{-4, 4}, {0, 2}}]}];
  sim = Table[{x,Exp[-(x - 2)^2]*RandomReal[{error1, error2}] + 
  Exp[-(x + 2)^2]*RandomReal[{error3, error4}]}, 
  {x,RandomReal[{-4, 4}, Points]}];
  G3 = ListPlot[sim, PlotStyle -> Thick];
  Show[If[bands, {G3, G2, G1}, G3], PlotRange -> {{-4, 4}, {0, 2.2}}],
  {{Points, 500}, 100, 1500, Appearance -> "Labeled"},
  {{error1, 1}, 0, 1, Appearance -> "Labeled"},
  {{error2, 1}, 1, 2, Appearance -> "Labeled"},
  {{error3, 0.5}, 0, 1, Appearance -> "Labeled"},
  {{error4, 1.5}, 1, 2, Appearance -> "Labeled"},
  Button["new sim", {Clear@sim, sim}, ImageSize -> 100],
  {{bands, False, Style["Show \"confidence\" bands?", Bold, Red, 
  FontSize -> 16]}, {True, False}}]

Answer 2

This question illustrates one of the main reasons to use Quantile regression -- heteroscedasticity. (Most of the classical time series analysis methods assume homoscedasticity.)

These documents/posts further exemplify time series fitting and reconstruction of CDFs at different time points:

• "Estimation of conditional density distributions",

• "Time series analysis with Quantile regression".

Here is code to do faithful fitting over the OP data with a B-spline basis. The function QuantileRegressionFit can be used to do the fit with user selected functions.

Import["https://raw.githubusercontent.com/antononcube/\
MathematicaForPrediction/master/QuantileRegression.m"]

qs = Range[0.1, 0.9, 0.2];
AbsoluteTiming[qFuncs = QuantileRegression[data, 20, qs];]

(* {0.54569, Null} *)

Show[{ListLinePlot[
   Table[{#, rq[#]} & /@ Sort[data[[All, 1]]], {rq, qFuncs}], 
   PlotTheme -> "Detailed", PlotLegends -> qs], 
  Graphics[{Gray, Point[data]}, PlotRange -> All]}, Frame -> True, 
 AspectRatio -> 1]