跳入整合()？

Question

我需要计算几个积分。我知道R不是做这件事的合适软件，但是由于我在R中做了其他的事情，而且积分只是一维的，我认为这可能是没有问题的。

无论如何，当增加积分的上限时，使用stats包中的函数integrate()会导致跳转。情节如下：

​

​

但是，如果通过构建和而不是积分，或者使用来自cubature包的函数cubature，结果如下所示：

​

​

为了能够复制这一点，下面是代码。我知道可能有一个更容易的例子，但这实际上是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中注释的第二段

当在无限间隔上进行积分时，要显式地这样做，而不是仅仅使用一个大的数字作为端点。这增加了得到正确答案的机会--在无限区间上积分有限的任何函数在该区间的大多数情况下都必须接近于零。

(“笔记”的其余部分也值得阅读。)

如果我把你的上限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中增长对象；首先分配整个向量，而不是重复地附加它们。)