หนึ่งแสดงช่วงข้อผิดพลาดในการวัดสำหรับฮิสโตแกรมอย่างไร

Question

ฉันมีปริมาณทางกายภาพแบบสุ่มซึ่งมีข้อผิดพลาดในการวัดที่เกี่ยวข้องด้วย มีวิธีที่ดีในการแสดงข้อผิดพลาดในการวัดบนฮิสโตแกรมโดยที่แกน x เป็นปริมาณความสนใจแบบสุ่มหรือไม่? มีวิธีอื่นในการมองเห็นทั้งการกระจายของปริมาณและข้อผิดพลาดในการวัดบนกราฟิกหนึ่งภาพหรือไม่

Answer 1

นี่อาจจะดูน่าเกลียดกว่าด้วยฮิสโตแกรม แต่ถ้าคุณมีข้อมูลเพียงพอสำหรับตัวอย่างบูตสแตรปที่จะทำการประมาณตัวอย่างต้นฉบับได้ดีคุณก็สามารถประมาณการกระจายตัวของตัวอย่างฮิสโตแกรมของคุณได้อย่างมีประสิทธิภาพและใช้ข้อมูลนั้นเพื่อให้ได้แถบความเชื่อมั่น

นี่คือตัวอย่างของ KDEs ข้อมูล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
))