Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

enter image description here

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

enter image description here

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

share|improve this question
It's assuming homoscedasticity – Searke May 20 '13 at 22:58
@Searke, I'm rather confused by the plot on the page you linked to. Naively I would expect that there would be tighter bands at the sides than in the middle of the plot as there seems to be less variance at the sides. Am I misunderstanding confidence bands or is it behaving in a strange way? – Jonathan Shock May 21 '13 at 8:33
Heteroscedasticity can be complicated: no longer does it suffice to stipulate that all residuals have the same variance, so now you need a model not only for the data values, but also an equally complicated one for the variances of their residuals. You're in the realm of highly specialized coding. If you know what you're doing, MMA is a great platform for writing that code; otherwise, look to specialized platforms like R and hope someone has contributed a model that might be appropriate for your situation. – whuber May 21 '13 at 15:23
So if I understand correctly the bottom line is that if I collect real data, I can use the prediction bands only if the variance of the residuals is constant, which I can evaluate by looking at fit["FitResiduals"]. – Ziofil May 21 '13 at 15:51

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.

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

enter image description here

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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