我们可以使用Base R来查找曲线下95％的面积吗？

Question

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

更具体地说，我想从mode（绿色虚线）移向尾部，然后在覆盖95％的曲线区域时停止。所需的是这个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