在 integrate() 中跳转？

Question

我需要计算一些积分。我知道 R 不是执行此操作的正确软件，但由于我在 R 中执行了所有其他操作并且积分只是一维的，所以我认为这可能没问题。

无论如何，当增加积分上限时，使用统计包中的函数 integrate() 会导致跳跃。这是情节：

但是，如果我通过求和而不是求积分，或者使用 cubature 包中的函数 adaptIntegrate()，结果如下所示：

为了能够复制它，这里是代码。我知道可能有一个更简单的例子，但这实际上是 1:1 我面临的情况。使用 integrate():

v_p = 11269
d_p = seq(0, v_p, 100)
tau = 0.35
interest = 0.05
time_t = 1
mu = 0
sigma_p = 0.1099313
nu_p = 0.0101552
ts = c()
bc = c()
for (i in 1:length(d_p)){
  ts  = c(ts,integrate(function(x) tau*(x-v_p-interest*d_p[i])*dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), v_p+interest*d_p[i], 10000000)$value)
  bc = c(bc,integrate(function(y) nu_p * y * dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), 0, d_p[i])$value)
}
shv = ts + bc
plot(d_p,shv, main = 290801)
abline(v = 9079, col="red")

完全相同的代码，只是使用 adaptIntegrate():

v_p = 11269
d_p = seq(0, v_p, 100)
tau = 0.35
interest = 0.05
time_t = 1
mu = 0
sigma_p = 0.1099313
nu_p = 0.0101552
ts = c()
bc = c()
for (i in 1:length(d_p)){
  ts  = c(ts,adaptIntegrate(function(x) tau*(x-v_p-interest*d_p[i])*dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), v_p+interest*d_p[i], 10000000)$integral)
  bc = c(bc,adaptIntegrate(function(y) nu_p * y * dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), 0, d_p[i])$integral)
}
shv = ts + bc
plot(d_p,shv, main = 290801)
abline(v = 9079, col="red")

有人知道为什么会这样吗？

Answer 1

您应该阅读 ?integrate 中 注意事项 的第二段：

When integrating over infinite intervals do so explicitly, rather than just using a large number as the endpoint. This increases the chance of a correct answer – any function whose integral over an infinite interval is finite must be near zero for most of that interval.

（注释的其余部分也值得一读。）

如果我将您的上限 10000000 更改为 Inf，我会得到一个合理的结果。（但请注意，这类问题通常 永远不会 完全解决……有时您需要针对您的特定问题调整间隔、起始点数等。）

用 Inf 替换的核心代码，为清楚起见略微重写：

f1 <- function(x) tau*(x-v_p-interest*d_p[i])*
        dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t))
f2 <- function(y) nu_p * y *
      dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t))
for (i in 1:length(d_p)){
  ts  = c(ts,integrate(f1, v_p+interest*d_p[i], Inf)$value)
  bc = c(bc,integrate(f2, 0, d_p[i])$value)
}

（顺便说一下，您还应该避免在 R 中增加对象；首先分配整个向量而不是重复附加到它们。）

在 integrate() 中跳转？

Jumps in integrate()?

r

integral