3
$\begingroup$

Here is the plot I want to make using Mathematica 8.0

enter image description here:

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.

$\endgroup$
  • $\begingroup$ To obtain the confidence bands you'll also need the estimated covariances of the parameter estimators. How did you obtain the estimates? $\endgroup$ – JimB May 22 '16 at 19:52
  • $\begingroup$ Actually I don't have those. Only the upper and lower limits plus the mean values are known from a model. $\endgroup$ – Benjamin May 22 '16 at 20:00
  • 3
    $\begingroup$ 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.) $\endgroup$ – JimB May 22 '16 at 20:11
  • $\begingroup$ 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, $\endgroup$ – Benjamin May 22 '16 at 20:14
  • 1
    $\begingroup$ 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.) $\endgroup$ – JimB May 22 '16 at 20:38
4
$\begingroup$

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.

enter image description here

enter image description here

|improve this answer|||||
$\endgroup$
  • $\begingroup$ This actually worked. Thanks Chris, $\endgroup$ – Benjamin May 22 '16 at 22:31
6
$\begingroup$

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

Correlation matrix

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

Data and fit

Your way...

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

Ops result

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

Confidence bands ignoring covariances

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Hi Jim. This is truly appreciated. $\endgroup$ – Benjamin May 23 '16 at 17:34
  • $\begingroup$ 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? $\endgroup$ – Benjamin May 23 '16 at 17:38
  • $\begingroup$ 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. $\endgroup$ – JimB May 23 '16 at 19:26
  • $\begingroup$ 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, $\endgroup$ – Benjamin May 23 '16 at 23:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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