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:


*

*Use bootstrapping to obtain 1000 samples of my data

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

*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!). 
 A: 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)


