我们能否使用基本的R语言来找到曲线下95%的面积？

Question

我们能否使用基本的R语言来找到曲线下95%的面积？

rfunctionoptimizationdistributionintegrate

8

使用基础R，我想知道是否可以确定下面标记为posterior的曲线下的95%区域？

更具体地说，我想从mode（绿色虚线）向两侧移动，然后在覆盖曲线面积的95％时停止。如下图所示，期望的是这个95％区域的x轴值限制。

     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％区间。

- rnorouzian

2个回答

1

这是一个利用梯形法则的解决方案。您会注意到@Remko提供的解决方案要好得多，但是该解决方案希望增加一些教育价值，因为它阐明了如何将复杂的问题简化为简单的几何、算术和基本编程结构，如for循环。

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

- Joseph Wood

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Remko Duursma · Accepted Answer

对称分布

虽然 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))
}

现在像这样使用上面的代码。我使用了一个非常不对称的函数（只是假设我的dist实际上是一些复杂的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在接近零的长序列上会出现问题。同样，对于您的情况，您需要调整域（除非有更好的想法！）。