对 R 中的积分进行积分答案

【问题标题】：Integrate over an integral in R对 R 中的积分进行积分
【发布时间】：2016-03-29 05:09:51
【问题描述】：

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

∫₀^H [π(t) ∫_t^HA(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

我不确定自己做错了什么。

【问题讨论】：

标签： r integral

【解决方案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 发现，将向量传递给 lower 和/或 upper 是允许的，但只考虑了第一个值。当在更宽的区间内积分时，正交方案选择离原点较远的第一个点，导致您观察到的不直观的下降。

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

【讨论】：

【解决方案2】：

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

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

【讨论】：