¿Cómo se muestra el rango de error de medición para un histograma?

Question

Tengo una cantidad física aleatoria que también tiene un error de medición asociado. ¿Existe una buena manera de mostrar el error de medición en un histograma donde el eje x es la cantidad aleatoria de interés? Alternativamente, ¿hay alguna otra forma de visualizar tanto la distribución de cantidades como el error de medición en un gráfico?

Answer 1

Esto puede ser más feo con un histograma, pero si tiene suficientes datos para que una muestra de arranque haga un buen trabajo de aproximación a la muestra original, entonces puede estimar efectivamente la distribución de muestreo de su histograma y usarla para obtener bandas de confianza.

Aquí hay un ejemplo con KDE. Los datos xse extraen iid de una distribución Gamma y se muestran como el diagrama de alfombra en la parte inferior. Si solo ajustamos un solo KDE, obtendríamos la línea negra gruesa. Pero podemos volver a muestrear xuna y otra vez y ajustar un KDE en cada muestra y trazar eso, lo cual se hace en rojo. Luego podemos tomar los cuantiles de 2,5 % y 97,5 % de las densidades remuestreadas para cada punto para tener una idea de la variación de la estimación puntual KDE. Esto es muy similar a muestrear una y otra vez la distribución posterior de una variable aleatoria y obtener bandas de confianza observando los cuantiles posteriores.

Aquí está el código para este ejemplo:

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)

---

Aquí hay un ejemplo con datos discretos: generé algunas observaciones iid $\text{Pois}(\lambda=8.54)$ y ajusté un histograma. Luego volví a muestrear los datos una y otra vez y calculé el histograma para cada nuevo muestreo utilizando los mismos contenedores que el original. Las barras de error provienen de los cuantiles de 2,5 % y 97,5 % de los histogramas resultantes.

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