R语言含乘积项被积函数的积分求解问题

问题根源 #

你的代码报错核心原因是：integrate函数会向量化调用被积函数（即一次传入多个x值），但你用循环赋值时，(n-i)/(x+n-i)会返回一个向量，却要赋值给term_11[i+1]这个单个元素，导致长度不匹配，进而计算出错误的0值。

正确实现方式 #

方案1：适配向量化输入的函数 #

修改函数，对每个传入的x值单独计算乘积项，避免长度不匹配问题：

my_func <- function(x) {
  n <- 45
  alpha_const <- 1
  beta_const <- 1
  
  # 用sapply遍历每个x值，单独计算对应结果
  sapply(x, function(x_val) {
    # 计算乘积项：遍历i从0到n，求所有(n-i)/(x_val+n-i)的乘积
    term_prod <- prod( (n - 0:n) / (x_val + n - 0:n) )
    # 计算被积函数值
    (x_val^(alpha_const - 1)) * exp(-beta_const * x_val) * term_prod
  })
}

# 求解积分
integrate(my_func, 0, Inf)

方案2：对数求和优化数值稳定性 #

当n=45时，乘积项可能因数值过小出现下溢（直接变成0），用对数求和代替直接乘积能提升计算精度：

my_func_stable <- function(x) {
  n <- 45
  alpha_const <- 1
  beta_const <- 1
  
  sapply(x, function(x_val) {
    # 将乘积转换为对数求和，再指数化还原
    log_terms <- log( (n - 0:n) / (x_val + n - 0:n) )
    term_prod <- exp(sum(log_terms))
    (x_val^(alpha_const - 1)) * exp(-beta_const * x_val) * term_prod
  })
}

# 求解积分
integrate(my_func_stable, 0, Inf)

关键说明 #

• integrate要求被积函数支持向量输入，用sapply逐个处理x值是最直接的适配方式。
• 对数转换能避免多次乘积导致的数值精度丢失，尤其当x较大时，效果更明显。

内容的提问来源于stack exchange，提问作者Arnoneel Sinha