How does one show measurement error range for a histogram? I have some random physical quantity that also has measurement error associated with it. Is there a good way to show the measurement error on a histogram where the x-axis is the random quantity of interest? Alternatively, is there some other way to visualize both the distribution of quantities and the measurement error on one graphic?
 A: This might be uglier with a histogram, but if you have enough data for a bootstrap sample to do a good job of approximating the original sample, then you can effectively estimate the sampling distribution of your histogram and use that to get confidence bands.
Here's an example with KDEs. The data x is drawn iid from a Gamma distribution and is shown as the rug plot at the bottom. If we just fit a single KDE we'd get the heavy black line. But we can resample from x over and over and fit a KDE on each sample and plot that, which is done in red. We can then take the 2.5% and 97.5% quantiles of the resampled densities for each point to get a sense of the variation of the point estimate KDE. This is very similar to sampling from the posterior distribution of a random variable over and over and getting confidence bands by looking at the posterior quantiles.

Here's the code for this example:
set.seed(1)
n <- 500
x <- rgamma(n, 2.34, 5.6)
d <- density(x)

nboot <- 5000
bootdat <- replicate(nboot, sample(x, n, TRUE))
dens <- apply(bootdat, 2, function(x) density(x)$y)
plot(0,0,col="white", xlim=range(d$x), ylim=c(0, max(d$y)*1.25), xlab="x", ylab="Density",
     main="Density estimate with bootstrap estimates")
apply(dens, 2, function(y) lines(y~d$x, col=rgb(red=1, green=0, blue=0, alpha=0.05)))
lines(d$y~d$x, lwd=3)  # the point estimate KDE

# computing and plotting the density quantiles
q <- apply(dens, 1, quantile, probs=c(.025, .975))
apply(q, 1, function(v) lines(v~d$x, col="blue", lwd=2, lty=2))
legend("topright", c("Point estimate", "Bootstrap estimate", "Bootstrap quantile"), col=c("black", "red", "blue"), bty="n", lty=c(1,1,2))
rug(x)


Here's an example with discrete data: I generated some iid $\text{Pois}(\lambda=8.54)$ observations and fit a histogram. I then resampled the data over and over and computed the histogram for each resample using the same bins as the original. The error bars come from the 2.5% and 97.5% quantiles of the resulting histograms.

set.seed(1)
sum_norm <- function(x) x / sum(x)
n <- 500
x <- rpois(n, 8.54)
h <- hist(x, 10, plot=FALSE)
h$counts <- sum_norm(h$counts)  # because `freq` ignored if `plot=FALSE`

nboot <- 5000
bootdat <- replicate(nboot, sample(x, n, TRUE))
hists <- apply(bootdat, 2, function(x) sum_norm(hist(x, breaks=h$breaks, plot=FALSE)$counts))

plot(h, ylim=range(hists), main = "Histogram with bootstrapped error bars", ylab = "Density")
q <- apply(hists, 1, quantile, probs=c(.025, .975))
midpts <- (h$breaks[-1] + h$breaks[-length(h$breaks)]) / 2
invisible(Map(
  function(y_lb, y_up, xpt)
    arrows(xpt, y_lb, xpt, y_up, col="red", code=3, angle=90, length=.05),
  q[1,], q[2,], midpts
))

