对 R 中的积分进行积分

Question

我想在 R 中解决以下问题：

∫₀^H [π(t) ∫_t^H A(x) dx ] dt

其中 π(t) 是先验，A(x) 是下面定义的 A 函数。

prior <- function(t) dbeta(t, 1, 24)
A     <- function(x) dbeta(x, 1, 4)
expected_loss <- function(H){
  integrand     <- function(t) prior(t) * integrate(A, lower = t, upper = H)$value
  loss          <- integrate(integrand, lower = 0, upper = H)$value
  return(loss)
}

由于 π(t), A(x) > 0，expected_loss(.5) 应小于 expected_loss(1)。但这不是我得到的：

> expected_loss(.5)
[1] 0.2380371
> expected_loss(1)
[1] 0.0625

我不确定我做错了什么。

Answer 1

在您的 integrand 中，lower = t 未矢量化，因此对集成的调用未达到您的预期*。矢量化 t 解决了这个问题，

expected_loss <- function(H){
  integrand <- function(t) prior(t) * integrate(A, lower = t, upper = H)$value
  vint <- Vectorize(integrand, "t")
  loss <- integrate(vint, lower = 0, upper = H)$value
  return(loss)
} 

expected_loss(.5)
# [1] 0.7946429
expected_loss(1)
# [1] 0.8571429

*：仔细观察 integrate 发现将向量传递到下部 and/or 上部是默默允许的，但只考虑第一个值。当在更宽的区间内积分时，正交方案选择了离原点更远的第一个点，导致您观察到的不直观的减少。

向 r-devel 报告此行为后，this user-error will now be caught by integrate 感谢 Martin Maechler (R-devel)。

Answer 2

在这种特殊情况下，您不需要 Vectorize，因为 dbeta 的积分已经通过 pbeta 在 R 中实现。试试这个：

prior <- function(t) dbeta(t, 1, 24)
#define the integral of the A function instead
Aint     <- function(x,H) pbeta(H, 1, 4) - pbeta(x,1,4)
expected_loss <- function(H){
  integrand<-function(x) Aint(x,H)*prior(x)
  loss          <- integrate(integrand, lower = 0, upper = H)$value
  return(loss)
}
expected_loss(.5)
#[1] 0.7946429
expected_loss(1)
#[1] 0.8571429

对 R 中的积分进行积分

Integrate over an integral in R

r

integral