3
$\begingroup$

NonlinearModelFit gave a bad fit to this data, but FindDistribution gave a good fit. This code does the fit and shows the results:

data= {0.0228278, 0.0235875, 0.0258227, 0.0281474, 0.0299132, 0.0300756, 
0.0301485, 0.0302263, 0.0306889, 0.030902, 0.0330661, 0.0357966, 
0.0361814, 0.0376612, 0.0402447, 0.0429684, 0.0442914, 0.0483349, 
0.0518819, 0.0529263, 0.0545642, 0.056154, 0.0591044, 0.0621301, 
0.0623543};
FindDistribution[data, 10, All, TargetFunctions -> "Continuous"]
Show[Histogram[data, Automatic, "PDF", Frame -> True], 
Plot[PDF[LogNormalDistribution[-3.266, 0.3239], x], {x, 0, 0.1}, 
PlotStyle -> Red, PlotRange -> All]]

Now I want to add 60% upper and lower confidence bands to the figure. How can this be done?

$\endgroup$
6
  • $\begingroup$ NonlinearModelFit is for regression and FindDistribution is for estimating distributions from a random sample. $\endgroup$
    – JimB
    Jan 17, 2018 at 0:09
  • $\begingroup$ There's an example in the Applications -> Confidence Intervals section of the PDF docs which might be useful. $\endgroup$ Jan 17, 2018 at 0:09
  • $\begingroup$ Would NonlinearModelFit with the same distribution as the fit function give the same result as FindDistribution? $\endgroup$ Jan 17, 2018 at 17:25
  • $\begingroup$ What you have is a random sample from some probability distribution and you want to estimate a probability density function. That is a totally different situation than performing a regression on data points that just happen to have the shape of a probability density function. Your choice of "60%" is also way out of the ordinary which suggests to me - and maybe wrongly - that you really ought to talk to a statistician. $\endgroup$
    – JimB
    Jan 17, 2018 at 18:38
  • $\begingroup$ I don't know any statisticians, but will try to find one. In the meantime, can you please elaborate on your comment? E.G., what is an ordinary confidence band? $\endgroup$ Jan 17, 2018 at 21:42

1 Answer 1

3
$\begingroup$

Perhaps, you can do bootstrap to get quick-and-dirty point estimates for standard deviation of PDF[dist, x]:

edist = EstimatedDistribution[data, LogNormalDistribution[μ, σ]];
bootstrap = PDF[EstimatedDistribution[#, LogNormalDistribution[μ, σ]], x] & /@ 
   RandomChoice[data, {100, IntegerPart[2 Length[data]/3]}];

plt = Plot[{Mean[bootstrap], Mean[bootstrap] - 1.96 StandardDeviation[bootstrap], 
    Mean[bootstrap] + 1.96 StandardDeviation[bootstrap]}, 
   {x, 0, 0.1}, Filling -> {1 -> {2}, 1 -> {3}}, FillingStyle -> Opacity[.5, Yellow]];

Show[Histogram[data, Automatic, "PDF", Frame -> True], 
 Plot[PDF[edist, x], {x, 0, 0.1}, PlotStyle -> Red], plt,  PlotRange -> All]

enter image description here

Change 1.96 to Quantile[NormalDistribution[], .6] to get

enter image description here

or, using bootstrap quantiles,

plt2 = Plot[{Median[bootstrap], Quantile[bootstrap, .40], Quantile[bootstrap, .60]}, 
  {x, 0, 0.1}, Filling -> {1 -> {2}, 1 -> {3}}, 
   FillingStyle -> Opacity[.5, Pink], PlotPoints -> 200];

Show[Histogram[data, Automatic, "PDF", Frame -> True], 
 Plot[PDF[edist, x], {x, 0, 0.1}, PlotStyle -> Red], plt2, PlotRange -> All]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks! I am unaware of "bootstrap," can you please give me some references? Also, does Quantile[NormalDistribution[], .6] give the 60% confidence bands? What does 1.96 give? $\endgroup$ Jan 17, 2018 at 17:23
  • $\begingroup$ @MichaelB.Heaney, wikipedia > Bootstrapping and Efron article. Re 1.96 it is Quantile[NormalDistribution[], .975] ( the "z-score" for two-sided 95% confidence interval.) $\endgroup$
    – kglr
    Jan 17, 2018 at 17:40

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.