Plot two curves and confidence bands representing the uncertainties in their parameters

Here is the plot I want to make using Mathematica 8.0

:

CDDHybrid[s_, a1_, b1_, c1_] := 10^(a1 + b1*c1 - b1*s)*Exp[-10^(s - c1)];
CDDHybridmean[s_] := 10^(-4.9684 - 0.88*s)*Exp[-10^(s - 20.82)];
CDD40[s_, a2_, b2_, c2_] := 10^(a2 + b2*c2 - b2*s)*Exp[-10^(s - c2)];
CDD40mean[s_] := 10^(-11.0765 - 0.57*s)*Exp[-10^(s - 20.55)];

With[{a1 = Interval[-23.29 + .04 {-1, 1}],
b1 = Interval[0.88 + .04 {-1, 1}],
c1 = Interval[20.82 + .01 {-1, 1}],
a2 = Interval[-22.79 + .06 {-1, 1}],
b2 = Interval[0.57 + .07 {-1, 1}],
c2 = Interval[20.55 + .03 {-1, 1}]},
LogPlot[{Min[CDDHybrid[s, a1, b1, c1]], CDDHybridmean[s],
Max[CDDHybrid[s, a1, b1, c1]], Min[CDD40[s, a2, b2, c2]],
CDD40mean[s], Max[CDD40[s, a2, b2, c2]]},
{s, 19, 23}, PlotRange -> {10^-28, 10^-20},
Filling -> {1 -> {3}, 4 -> {6}},
FillingStyle -> {LightBlue, LightRed},
PlotStyle -> {Blue, Blue, Blue, Red, Red, Red, Brown, Green},
Frame -> True,
FrameLabel -> {Style["X [some unit]", FontSize -> 26], Style["Y [some unit]", FontSize -> 26]},
FrameTicksStyle -> Directive[FontSize -> 26]]]


I am getting the plots correctly but it seems that the code doesn't seem to recognize that two fillings have different colors. I don't know how to fix the issue. Attached is my plot for reference.

• To obtain the confidence bands you'll also need the estimated covariances of the parameter estimators. How did you obtain the estimates? – JimB May 22 '16 at 19:52
• Actually I don't have those. Only the upper and lower limits plus the mean values are known from a model. – Benjamin May 22 '16 at 20:00
• If you don't have the data either and you only have 1 digit for the standard errors (I assume that the plus-or-minus numbers are standard errors), then you shouldn't produce intervals or confidence bands which can only be misleading and maybe dangerous if you're developing weapons systems or similar things. (If your boss or major professor is telling you to do this, then you have my sympathy.) – JimB May 22 '16 at 20:11
• Could you please let me know if there were data with their error bars (instead of a functional form), what would be the syntax then? I might need this in future too. Thanks, – Benjamin May 22 '16 at 20:14
• I'm sure your data don't have error bars. You have a (schizophrenic) function (i.e., one that uses both 10^x and Exp) that likely is fit after taking logs and there might be some experimental design such as obtaining CDDHybrid and CDD40 on the same experimental unit. So the answer is "It depends." It depends on what the error structure is: does the variability about the curve stay constant over all values of s? Does it increase with increasing s? Is the fit good? (No point in having confidence bands for a poor fit.) – JimB May 22 '16 at 20:38

Could use

Filling -> {4 -> {{1}, LightRed}, 1 -> {{3}, LightBlue}, 3 -> {{6}, LightRed}}


and omit FillingStyle.

Or

Filling -> {1 -> {3}, 4 -> {6}},
FillingStyle -> {{LightBlue, Opacity[0.5]}, {LightRed, Opacity[0.5]}}


which respectively produce the plots below.

• This actually worked. Thanks Chris, – Benjamin May 22 '16 at 22:31

Below is an extended comment showing the effects that can happen when one ignores the correlations among parameter estimators (either on purpose or because journal authors or editors don't know enough about statistics to include necessary information).

First I create some data with the same curve form as the function CDDHybrid and a set of predictor values and error about the curve that results in very similar values of the reported standard errors. Then the confidence bands are created.

(* Function to be considered *)
CDDHybrid[s_, a_, b_, c_] := 10^(a + b*c - b*s)*Exp[-10^(s - c)];

(* Generate some data *)
n = 20; (* Sample size *)
(* Predictor variables *)
x = Table[19. + 2.8 (i - 1)/(n - 1), {i, n}];
(* Responses *)
SeedRandom[1234]
y = Flatten[
Table[Log[CDDHybrid[x[[i]], -23.29, 0.88, 20.82]] +
RandomVariate[NormalDistribution[0, 0.2], 1], {i, n}]];

(* Find estimates *)
nlm = NonlinearModelFit[Transpose[{x, y}], (a + b*c - b*s) Log[10] - 10^(s - c),
{{a, -23.29}, {b, 0.88}, {c, 20.82}}, s];
estimates = {a, b, c} /. nlm["BestFitParameters"]
(* {-23.2725, 0.866555, 20.8164} *)
corr = nlm["CorrelationMatrix"] // MatrixForm


Note the very high correlations. And below we see that the estimated standard errors are similar to the reported values.

(* Standard errors of parameter estimators *)
vcov = nlm["CovarianceMatrix"];
se = Diagonal[vcov]^0.5
(* {0.04771502502282548, 0.03662121809566094 ,0.01209891148269071} *)


Now plot the confidence bands...

confidenceBands = nlm["MeanPredictionBands"];
Show[LogPlot[Exp[confidenceBands], {s, 19, 23},
PlotRange -> {{19, 23}, {10^-28, 10^-20}},
Filling -> {3 -> {1}}, FillingStyle -> {LightBlue},
PlotStyle -> {Blue, Blue, Blue},
Frame -> True, FrameLabel -> {Style["X [some unit]", FontSize -> 26],
Style["Y [some unit]", FontSize -> 26]},
FrameTicksStyle -> Directive[FontSize -> 26]],
ListLogPlot[Transpose[{x, Exp[y]}],
PlotStyle -> {Red, PointSize[0.03]}]]


CDDHybridmean[s_] := 10^(estimates[[1]] + estimates[[2]]*estimates[[3]] -
estimates[[2]]*s)*Exp[-10^(s - estimates[[3]])]
With[{a1 = Interval[estimates[[1]] + se[[1]] {-1, 1}],
b1 = Interval[estimates[[2]] + se[[2]] {-1, 1}],
c1 = Interval[estimates[[3]] + se[[3]] {-1, 1}]},
LogPlot[{Min[CDDHybrid[s, a1, b1, c1]], CDDHybridmean[s],
Max[CDDHybrid[s, a1, b1, c1]]}, {s, 19, 23},
PlotRange -> {10^-28, 10^-20}, Filling -> {1 -> {3}},
FillingStyle -> {LightBlue}, PlotStyle -> {Blue, Blue, Blue},
Frame -> True,
FrameLabel -> {Style["X [some unit]", FontSize -> 26],
Style["Y [some unit]", FontSize -> 26]},
FrameTicksStyle -> Directive[FontSize -> 26]]]


Update in response to a comment...

If one is going to assume (which should be "able to assume" rather than just "willing to assume") that the covariances among the parameter estimators are zero, then the Delta Method (or Propagation of Error Method) can be used under appropriate conditions (which I won't go into) to estimate variances and confidence bands. Assuming covariances are zero does not necessarily make the confidence bands conservative.

(* Function being fit *)
f = (a + b (c - s)) Log[10] - 10^(s - c);

(* Estimates of parameters and associated variances *)
variances = {0.04, 0.04, 0.01}^2;
estimates = {a -> -23.29 , b -> 0.88, c -> 20.82};

(* Use Delta/Propagation of error method to approximate variance of mean at a value of s *)
(* http://en.wikipedia.org/wiki/Propagation_of_uncertainty *)
(* http://en.wikipedia.org/wiki/Delta_method *)
(* (But unwisely assuming covariances are all zero) *)
d = D[f, {{a, b, c}}];
se = ((d^2). variances /. estimates)^0.5
mean = f /. estimates

(* Plot results *)
LogPlot[{Exp[mean - 1.96 se], Exp[mean], Exp[mean + 1.96 se]}, {s, 19, 23},
PlotRange -> {10^-28, 10^-20}, Filling -> {1 -> {3}},
FillingStyle -> {LightBlue},
PlotStyle -> {Blue, Blue, Blue}, Frame -> True,
FrameLabel -> {Style["X [some unit]", FontSize -> 26],
Style["Y [some unit]", FontSize -> 26]},
FrameTicksStyle -> Directive[FontSize -> 26]]


• Hi Jim. This is truly appreciated. – Benjamin May 23 '16 at 17:34
• I have to contact the authors of the paper and get their actual data and errors to predict better values for confidence bands. I am not going to trust my big bands then. Can we say that if off-diagonal elements are negligible then I will be getting such a wide band? – Benjamin May 23 '16 at 17:38
• I'v updated the extended comment to show the results of applying the Delta Method but assuming that the covariances are zero. But I have to ask: With an estimate of $a$ being -23.29 and the standard error being only 0.04, isn't that awfully precise? (I don't know the subject matter or how many samples were take so I have to ask.) And, no, it is your approach that makes things so wide. That's not how it should be done. – JimB May 23 '16 at 19:26
• Yes. The results is really precise. I have to admit that. But since I had no access to the actual data except their best fit curve, I thought if I proceed with this method (in Mathematica) I should get small scatter which turned out to be not the case and proved that the method is incorrect to what I expected. Again thanks for being so concerned, – Benjamin May 23 '16 at 23:47