Como mostrar a faixa de erro de medição para um histograma?

Question

Eu tenho uma quantidade física aleatória que também possui um erro de medição associado a ela. Existe uma boa maneira de mostrar o erro de medição em um histograma onde o eixo x é a quantidade aleatória de interesse? Alternativamente, existe alguma outra maneira de visualizar a distribuição de quantidades e o erro de medição em um gráfico?

Answer 1

Isso pode ser mais feio com um histograma, mas se você tiver dados suficientes para uma amostra bootstrap fazer um bom trabalho de aproximação da amostra original, poderá estimar efetivamente a distribuição de amostragem de seu histograma e usá-lo para obter faixas de confiança.

Aqui está um exemplo com KDEs. Os dados xsão desenhados a partir de uma distribuição Gama e são mostrados como o gráfico de tapete na parte inferior. Se apenas encaixarmos um único KDE, obteremos a linha preta pesada. Mas podemos reamostrar xrepetidamente e ajustar um KDE em cada amostra e plotar isso, o que é feito em vermelho. Podemos então obter os quantis de 2,5% e 97,5% das densidades reamostradas para cada ponto para obter uma noção da variação da estimativa pontual KDE. Isso é muito semelhante à amostragem da distribuição posterior de uma variável aleatória repetidas vezes e à obtenção de faixas de confiança observando os quantis posteriores.

Aqui está o código para este exemplo:

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)

---

Aqui está um exemplo com dados discretos: Gerei algumas observações iid $\text{Pois}(\lambda=8.54)$ e ajustei um histograma. Em seguida, fiz uma nova amostra dos dados várias vezes e calculei o histograma para cada nova amostra usando as mesmas caixas do original. As barras de erro vêm dos quantis de 2,5% e 97,5% dos 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
))