ヒストグラムの測定誤差範囲をどのように表示しますか?

Aug 19 2020

測定誤差も伴うランダムな物理量があります。x軸が対象のランダムな量であるヒストグラムに測定誤差を表示する良い方法はありますか?あるいは、1つのグラフィックで量の分布と測定誤差の両方を視覚化する他の方法はありますか?

回答

2 jld Aug 19 2020 at 12:06

これはヒストグラムでは見苦しいかもしれませんが、ブートストラップサンプルが元のサンプルを適切に近似するのに十分なデータがある場合は、ヒストグラムのサンプリング分布を効果的に推定し、それを使用して信頼帯を取得できます。

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