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

Question

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

更具体地说，我想从 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)

P.S换句话说，如果这样的间隔是最短的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 在处理接近零的长序列时会遇到问题。同样，对于您的情况，您将不得不调整域（除非有人有更好的主意！）。

Answer 2

这是一个利用 Trapezoidal rule 的解决方案。您会注意到@Remko 提供的解决方案要好得多，但是这个解决方案有望增加一些教学价值，因为它阐明了如何将复杂的问题简化为简单的几何、算术和基本编程结构，例如 for loops.

findXVals <- function(lim, p) {
    ## (1/p) is the precision

    ## area of a trapezoid
    trapez <- function(h1, h2, w) {(h1 + h2) * w / 2}

    yVals <- posterior((1:(p - 1))/p)
    m <- which.max(yVals)
    nZ <- which(yVals > 1/p)

    b <- m + 1
    e <- m - 1
    a <- f <- m

    area <- 0
    myRng <- 1:(length(nZ)-1)
    totArea <- sum(trapez(yVals[nZ[myRng]], yVals[nZ[myRng+1]], 1/p))
    targetArea <- totArea * lim

    while (area < targetArea) {
        area <- area + trapez(yVals[a], yVals[b], 1/p) + trapez(yVals[e], yVals[f], 1/p)
        a <- b
        b <- b + 1
        f <- e
        e <- e - 1
    }

    c((a - 1)/p, (f + 1)/p)
}

findXVals(.95, 10^5)
[1] 0.66375 0.48975

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

Can we use Base R to find the 95% of the area under a curve?

r

function

optimization

distribution

integrate

对称分布

分布不对称