# Adding uncertainty range to probability density function using bootstrapping

I'm hoping this isn't a ridiculous question! Here goes...

I would like to visualise the range of uncertainty in a probability density function fit to observed data using Maximum Likelihood estimation. I have included a sketch to illustrate what I am aiming for (note this is not based on real data, just a sketch). The steps I have used so far are:

1. Use bootstrapping to obtain 1000 samples of my data

2. Use maximum likelihood estimation to fit a three paramter inverse gamma probability density function (aka Pearson V distribution) to each sample from step 1 .

3. Taking the median of each parameter value from step 2 to plot the "median" distribution. Ditto for the min, max etc.

However, I think I may be going wrong at step 3, as the median of one parameter value does not necessarily correspond to the median of the other two. Is there a more statistically correct way to go about this?

Thanks in advance for any help - Cross Validated has become my go-to for so many questions over the last year (I just wish I could answer more rather than just ask!).

• Normally, people don't use bands around PDFs as in your illustration, because they are unrealistic: the upper band would have total probability greater than $1$ while the lower would have total probability less than $1$. A standard technique draws bands around the cumulative distribution functions. Would such a technique be acceptable to you? – whuber Jul 15 '13 at 14:32
• Thanks @Whuber, that makes sense to me that the probabilities would not sum to 1. It will take me a little time to take in the Probability Box concept but I think this could be helpful. Just thinking aloud, I suppose one way of visualising the range of pdfs would be to plot them all in one figure. Thanks again! – Faith Jul 15 '13 at 14:41

I'll go out on a limb and disagree with @whuber here. I don't think there's anything wrong with putting bands around pdfs, as long as you understand what they are: pointwise errors. It's like the similar confidence bands around smooth curves for additive regression models. If the bands are wide enough, they make it look like the curve could be far away from the center of the data, but people generally understand what's going on.

What I'd do is obtain the confidence bands and median of your density function directly from the bootstrap samples. You have 1000 estimated functions, so you can obtain your quantiles as the 2.5%, 50% and 97.5% points of the distribution of these functions at any given x.

The main downside is that there's no guarantee the resulting quantiles will actually be members of your parametric density family. Depending on your application, this may or may not be a problem. I doubt it would be, since the bootstrap is generally used only to get estimates of error, and you already have a parametric result from your original sample.

Here's a simple example using the gamma density and moment estimators. You can easily substitute the inverse gamma and maximum likelihood.

# the data
x <- rgamma(100, 2, scale=1)

# get 1000 bootstrapped estimates of the density, over a regular grid
xs <- seq(0, 8, len=201)
ests <- sapply(1:1000, function(i) {
xi <- sample(x, size=length(x), replace=TRUE)
vi <- var(xi)
mi <- mean(xi)
dgamma(xs, shape=mi^2/vi, scale=vi/mi)
})

# plot the individual estimates
plot(xs, ests[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(ests))
for(i in 2:ncol(ests)) lines(xs, ests[, i], col=rgb(.6, .6, .6, .1))

# get the 2.5%, 50% and 97.5% quantiles of the density estimate at each grid point
quants <- apply(ests, 1, quantile, c(0.025, 0.5, 0.975))
lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[2, ], col="darkred", lwd=2) • +1. And by substituting dgamma by pgamma, the very same script can be used to draw bootstrap bands around the cumulative distribution. – COOLSerdash Jul 15 '13 at 18:24
• +1 - I might also suggest since the OP is fitting ML estimates (and not say just estimating a KDE on the bootstrapped samples) an additional plot might be scatterplots of the bootstrapped parameter estimates. If a reviewer had a problem with this an easy fix is to just remove the quantile lines, the individual bootstrap estimates greyed out are valid PDF's and provide a good visual cue of the uncertainty. – Andy W Jul 15 '13 at 18:37
• Thanks hong-ooi for providing such clear code, this works really well. I mainly want this to give a fairly intuitive visualisation of the uncertainty of the pdf fit to the data (as I have several different data sets, each with a different fit this is particularly good for illustrating how certain data sets are more or less variable). I think I just need to be clear that these aren't actually "error" bars, as whuber states. Thanks also @Andy-W, I'm producing boxplots of the parameter estimates. – Faith Jul 17 '13 at 9:27