我们可以使用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)

PS换句话说,非常希望这样的间隔是可能的最短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建议找到最短的间隔.

解决方案中重要的是a的定义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长序列的近零有问题.同样,对于您的情况,您将不得不调整域名(除非有人有更好的想法!).