【问题标题】:R integrate: returns wrong solution (is using wrong quadrature points?)R 积分:返回错误的解(是否使用了错误的交点?)
【发布时间】:2015-02-02 12:40:37
【问题描述】:

我在 R 中有一个函数,我正在尝试集成它,但对于函数参数的某些(极端)值,integrate 返回不正确的解决方案。我相信问题可能是integrate 为其中一些极值选择了不正确的正交点,但首先我将提供演示该问题。

我希望集成的功能如下。

integrandFunc_F <- function(x, func_u, func_u_lowerBar, 
  func_u_upperBar, func_mean_v, func_sigma_v, func_sigma_epsilon, 
  func_sigma_y, func_gamma, func_rho) {
#print(x);
p <- 1 - pnorm(func_u_upperBar,x,func_sigma_y);
q <- pnorm(func_u_lowerBar,x,func_sigma_y);
p <- p*(1-func_rho); q <- q*(1-func_rho);
alpha <- ifelse(func_gamma*(p+q) == 0, 0, pmax((func_gamma*p-q)/(func_gamma*(p+q)), 0));
g <- ifelse(x > func_u, dnorm(x,func_mean_v,sqrt(func_sigma_v^2 + func_sigma_epsilon^2))/(1-pnorm(func_u,func_mean_v,sqrt(func_sigma_v^2 + func_sigma_epsilon^2))), 0);
output <- alpha*g;
output
}

当我尝试计算以下内容时,我得到了 1 的正确解:

integrate(integrandFunc_F, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
  func_u_upperBar= 8, func_mean_v= 30, func_sigma_v= .1, func_sigma_epsilon= 2, 
  func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)

但是,当我尝试计算以下内容时,我得到了 0 的错误解:

integrate(integrandFunc_F, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
  func_u_upperBar= 8, func_mean_v= 50, func_sigma_v= .1, func_sigma_epsilon= 2, 
  func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)

上面我表示我相信这个问题可能与正交点的选择有关。如果您在上面的函数中取消注释#print(x),您可以看到在func_mean_v = 30 的情况下,integrate 定位在相对较大/接近 30 的正交点上。但是,在func_mean_v=50 的情况下,经过几次迭代 @ 987654331@ 选择接近 0 的正交点。接近 0 的正交点不适合评估此函数,该函数包含均值在func_mean_v. 的正态分布

关于如何解决这个问题的任何想法?为什么integrate 在某些情况下会迭代到接近 0 的正交点?请注意,func_mean_v = 30func_mean_v = 50 的选择无疑是这个函数的极端参数,但是我需要能够正确计算这种情况。

【问题讨论】:

  • 您是否尝试过降低integrate收敛的容差?没有适用于所有被积函数的自适应数值积分的一般规则(它怎么知道被积函数在 -Inf 和 +Inf 之间的非零位置?)所以你有时need to help it
  • 我尝试降低rel.tol 并增加subdivisions 的数量。两者都没有解决问题。有没有办法为integrate 提供一组正交点?对于正态分布,不难知道被积函数在哪里非零。
  • 如果您知道被积函数在哪里非零,则积分变量的移位和可选重新缩放是imho的方法;这将极大地帮助integrate
  • integrate() 不能使用您自己的正交点:它使用特定的(高斯)正交规则,因此权重仅适用于相应的节点。

标签: r integrate numerical-integration


【解决方案1】:

您可以将积分变量移动到峰值的中心,

wrapper <- function(x, func_mean_v, ...)
   integrandFunc_F(x+func_mean_v, func_mean_v=func_mean_v, ...)


integrate(wrapper, rel.tol = 1e-8, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
          func_u_upperBar= 8, func_mean_v= 50, func_sigma_v= .1, func_sigma_epsilon= 2, 
          func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)
# 1 with absolute error < 1.3e-09

【讨论】:

  • 谢谢@baptiste。您可以编辑以显示为您提供解决方案 1 的确切代码吗?包装器是在其中定义了 integrandFunc_F 的函数吗?还是分别定义?
  • 准确的代码!您可能会对包装器缺少 {} 感到困惑,它们对于单行代码不是必需的。
  • 谢谢!所以我保留了我之前对integrandFunc_F 的原始定义,然后像您在此处所做的那样定义了函数wrapper
猜你喜欢
  • 2021-07-23
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2018-01-24
  • 2014-09-26
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多