Bagaimana seseorang menunjukkan rentang kesalahan pengukuran untuk histogram?

Aug 19 2020

Saya memiliki beberapa kuantitas fisik acak yang juga memiliki kesalahan pengukuran yang terkait dengannya. Adakah cara yang baik untuk menunjukkan kesalahan pengukuran pada histogram di mana sumbu x adalah kuantitas acak yang diinginkan? Atau, adakah cara lain untuk memvisualisasikan distribusi kuantitas dan kesalahan pengukuran pada satu grafik?

Jawaban

2 jld Aug 19 2020 at 12:06

Ini mungkin lebih buruk dengan histogram, tetapi jika Anda memiliki cukup data untuk sampel bootstrap untuk melakukan pekerjaan yang baik dalam mendekati sampel asli, maka Anda dapat memperkirakan distribusi sampling histogram Anda secara efektif dan menggunakannya untuk mendapatkan pita kepercayaan.

Berikut contoh KDE. Data xdiambil dari distribusi Gamma dan ditampilkan sebagai plot karpet di bagian bawah. Jika kita hanya memasukkan satu KDE kita akan mendapatkan garis hitam pekat. Tetapi kita dapat xmengambil sampel berulang-ulang dan memasukkan KDE pada setiap sampel dan plot itu, yang dilakukan dengan warna merah. Kami kemudian dapat mengambil kuantitas 2,5% dan 97,5% dari kepadatan sampel ulang untuk setiap titik untuk mendapatkan gambaran variasi perkiraan titik KDE. Ini sangat mirip dengan pengambilan sampel dari distribusi posterior variabel acak berulang kali dan mendapatkan pita kepercayaan dengan melihat kuantil posterior.

Berikut kode untuk contoh ini:

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)

Berikut adalah contoh dengan data diskrit: Saya membuat beberapa iid $ \ text {Pois} (\ lambda = 8.54) $ observasi dan menyesuaikan histogram. Saya kemudian mengambil sampel ulang data berulang-ulang dan menghitung histogram untuk setiap sampel ulang menggunakan tempat sampah yang sama seperti aslinya. Bar kesalahan berasal dari jumlah 2,5% dan 97,5% dari histogram yang dihasilkan.

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