gpt4 book ai didi

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

转载 作者:行者123 更新时间:2023-12-03 15:40:44 25 4
gpt4 key购买 nike

使用 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% 间隔,则是非常可取的。

enter image description here

最佳答案

对称分布

尽管 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)



enter image description here

不对称分布

在非对称分布的情况下,我们必须搜索满足 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对长序列的接近零有问题。同样,对于您的情况,您将不得不调整域(除非有人有更好的主意!)。

enter image description here

关于r - 我们可以使用 Base R 来找到曲线下 95% 的面积吗?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/45702886/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com