Как показать диапазон ошибок измерения на гистограмме?

Question

У меня есть некоторая случайная физическая величина, с которой также связана ошибка измерения. Есть ли хороший способ показать ошибку измерения на гистограмме, где по оси абсцисс отложена интересующая случайная величина? В качестве альтернативы, есть ли другой способ визуализировать как распределение величин, так и ошибку измерения на одном графике?

Answer 1

Это может быть уродливее с гистограммой, но если у вас достаточно данных для выборки начальной загрузки, чтобы хорошо аппроксимировать исходную выборку, то вы можете эффективно оценить распределение выборки вашей гистограммы и использовать это для получения доверительных интервалов.

Вот пример с KDE. Данные xвзяты из гамма-распределения и показаны в виде коврового графика внизу. Если мы просто поместим один KDE, мы получим жирную черную линию. Но мы можем повторно семплировать xснова и снова, подогнать KDE к каждому семплу и нанести его на график, что делается красным цветом. Затем мы можем взять 2,5% и 97,5% квантилей повторной выборки плотностей для каждой точки, чтобы получить представление об изменении точечной оценки KDE. Это очень похоже на многократную выборку из апостериорного распределения случайной величины и получение доверительных интервалов путем просмотра апостериорных квантилей.

Вот код этого примера:

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)

---

Вот пример с дискретными данными: я создал несколько наблюдений iid $ \ text {Pois} (\ lambda = 8.54) $ и подобрал гистограмму. Затем я повторно семплировал данные снова и снова и вычислял гистограмму для каждой повторной выборки, используя те же интервалы, что и оригинал. Планки погрешностей взяты из квантилей 2,5% и 97,5% полученных гистограмм.

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