在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，预期损失(.5)应该小于预期损失(1)。但这并不是我得到的结果:

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

我不确定我做错了什么。

- LJB

Answer 1

在您的被积函数中，lower = 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函数后发现可以将向量传递给上限和/或下限，但仅考虑第一个值。当在更广的区间上进行积分时，求积方案选择离原点更远的第一个点，导致你所观察到的不直观的减少。

在向r-devel报告了这种行为之后，由于Martin Maechler(R-devel)的贡献，此用户错误现在将被integrate函数捕获。

- baptiste

Answer 2

在这种情况下，您不需要进行向量化处理，因为在R中，通过使用pbeta，已经实现了对dbeta的积分。请尝试以下内容：

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

- nicola