r - 我们可以使用 Base R 来找到曲线下 95% 的面积吗？

Question

使用 Base R，我想知道是否可以确定曲线下 95% 的面积，表示为 posterior以下？

更具体地说，我想从 mode 移动(绿色虚线)朝向尾部，然后在我覆盖 95% 的曲线区域时停止。想要的 x 轴值是这个 95% 区域的限制，如下图所示？

     prior = function(x) dbeta(x, 15.566, 7.051) 
likelihood = function(x) dbinom(55, 100, x)
 posterior = function(x) prior(x)*likelihood(x)

mode = optimize(posterior, interval = c(0, 1), maximum = TRUE, tol = 1e-12)[[1]]

curve(posterior, n = 1e4)

附注 换句话说，如果这样的间隔是可能的最短 95% 间隔，则是非常可取的。

Answer 1

对称分布

尽管 OP 的示例并不完全对称，但它已经足够接近了 - 从那里开始很有用，因为解决方案要简单得多。

您可以使用 integrate 的组合和 optimize .我将此作为自定义函数编写，但请注意，如果您在其他情况下使用它，则可能需要重新考虑搜索分位数的界限。

# For a distribution with a single peak, find the symmetric!
# interval that contains probs probability. Search over 'range'.
f_quan <- function(fun, probs, range=c(0,1)){

  mode <- optimize(fun, interval = range, maximum = TRUE, tol = 1e-12)[[1]]

  total_area <- integrate(fun, range[1], range[2])[[1]]

  O <- function(d){
    parea <- integrate(fun, mode-d, mode+d)[[1]] / total_area
    (probs - parea)^2
  }
  # Bounds for searching may need some adjustment depending on the problem!
  o <- optimize(O, c(0,range[2]/2 - 1E-02))[[1]]

return(c(mode-o, mode+o))
}

像这样使用它，

f <- f_quan(posterior, 0.95)
curve(posterior, n = 1e4)
abline(v=f, col="blue", lwd=2, lty=3)

给



不对称分布

在非对称分布的情况下，我们必须搜索满足 P(a < x < b) = Prob 条件的两个点，其中 Prob 是某个期望的概率。由于有无限多个区间 (a,b) 满足这一点，因此 OP 建议找到最短的区间。

解决方案中重要的是定义 domain ，我们要搜索的区域(我们不能使用 -Inf, Inf ，因此用户必须将其设置为合理的值)。

# consider interval (a,b) on the x-axis
# integrate our function, normalize to total area, to 
# get the total probability in the interval
prob_ab <- function(fun, a, b, domain){
  totarea <- integrate(fun, domain[1], domain[2])[[1]]
  integrate(fun, a, b)[[1]] / totarea
}

# now given a and the probability, invert to find b
invert_prob_ab <- function(fun, a, prob, domain){

  O <- function(b, fun, a, prob){
    (prob_ab(fun, a, b, domain=domain) - prob)^2
  }

  b <- optimize(O, c(a, domain[2]), a = a, fun=fun, prob=prob)$minimum

return(b)
}

# now find the shortest interval by varying a
# Simplification: don't search past the mode, otherwise getting close
# to the right-hand side of domain will give serious trouble!
prob_int_shortest <- function(fun, prob, domain){

  mode <- optimize(fun, interval = domain, maximum = TRUE, tol = 1e-12)[[1]]

  # objective function to be minimized: the width of the interval
  O <- function(a, fun, prob, domain){
    b <- invert_prob_ab(fun, a, prob, domain)

    b - a
  }

  # shortest interval that meets criterium
  abest <- optimize(O, c(0,mode), fun=fun, prob=prob, domain=domain)$minimum

  # now return the interval
  b <- invert_prob_ab(fun, abest, prob, domain)

return(c(abest,b))
}

现在像这样使用上面的代码。我使用了一个非常不对称的函数(假设 mydist 实际上是一些复杂的 pdf，而不是 dgamma)。

mydist <- function(x)dgamma(x, shape=2)
curve(mydist(x), from=0,  to=10)
abline(v=prob_int_shortest(mydist, 0.9, c(0,10)), lty=3, col="blue", lwd=2)

在这个例子中，我将域设置为 (0,10)，因为显然间隔必须在某个地方。请注意，使用非常大的值(如 (0, 1E05))不起作用，因为 integrate对长序列的接近零有问题。同样，对于您的情况，您将不得不调整域(除非有人有更好的主意!)。