二参数贝叶斯可信区间

Question

如何确定多参数模型后验估计的 confidence/credibility 区间？

我可以分别得到每个参数的置信区间。（目前正在使用 bayestestR，但我不介意使用其他东西）

library(dplyr)
library(ggplot2)
library(bayestestR)

N <- 10000
#Posterior samples (random example)
p1 <- rnorm(N)
p2 <- p1 + rnorm(N)

df_post <- tibble(p1,p2)

describe_posterior(
  df_post,
  centrality = "median",
  test = c("p_direction", "p_significance")
)

##Yields:
# Summary of Posterior Distribution
# 
# Parameter |    Median |        95% CI |     pd |   ps
# -----------------------------------------------------
#   p1      |      0.02 | [-1.85, 2.08] | 50.64% | 0.46
#   p2      | -2.24e-03 | [-2.82, 2.78] | 50.04% | 0.47

我可以生成一个包含点和 2D 等高线的图，这为我提供了后验分布的视觉指示（尽管我不知道每个等高线代表什么百分比）：

ggplot(df_post, aes(x=p1, y=p2)) +
  geom_density_2d(size=1) +
  geom_point(size=0.1)

我的问题是，我如何实际计算（and/or 绘图）二维 X% 可信区间？

Answer 1

这是一个 base-R 绘图解决方案，它根据二维核密度估计绘制了 95% 的最高后验密度区域：

library(emdbook)
library(coda)
HPDregionplot(as.mcmc(df_post))
with(df_post, points(p1, p2, col = adjustcolor("black", alpha.f = 0.2)))

更多 smoothly/within ggplot:

library(ggplot2); theme_set(theme_bw())
## see function definition below
L <- with(df_post, get_hpd2d_levels(p1, p2))
gg0 <- ggplot(df_post, aes(p1, p2)) + geom_point(alpha=0.1) +
  geom_density_2d(breaks = L, colour="red")
print(gg0)

后验密度最高的区域是经典贝叶斯方法。如果您想深入探索，可以看看一些参数较少的方法来计算中心集（包图、功能深度等）。这类似于最高后验密度区域和基于分位数的可信区间之间的差异。

##' Get 2D highest posterior density levels corresponding to probability regions
##' @param x x-coordinate of samples
##' @param y y-coordinate
##' @param probs vector of probability levels
##' @param ... arguments for MASS::kde2d
##' @examples
##' dd <- data.frame(x=rnorm(1000), y=rnorm(1000))
##' get_hpd2d_levels(dd$x,dd$y)
##' gg2 <- ggplot(dd, aes(x,y)) + geom_density_2d(breaks=get_hpd2d_levels(dd$x,dd$y), colour="red")
##' print(gg2)
get_hpd2d_levels <- function(x, y, prob=c(0.9,0.95), ...) {
  post1 <- MASS::kde2d(x, y)
  dx <- diff(post1$x[1:2])
  dy <- diff(post1$y[1:2])
  sz <- sort(post1$z)
  c1 <- cumsum(sz) * dx * dy
  ## remove duplicates
  ## dups <- duplicated(sz)
  ## sz <- sz[!dups]
  ## c1 <- c1[!dups]
  levels <- sapply(prob, function(x) {
    approx(c1, sz, xout = 1 - x, ties = mean)$y
  })
  return(levels)
}

Answer 2

works well, but here's another suggestion (based on ) 也有效：

car::dataEllipse(p1, p2, levels=0.95)

二参数贝叶斯可信区间

Two-parameter Bayesian credibility Interval

r

bayesian