使用integrate与uniroot求解含Weibull、指数分布的二重积分及beta值

错误原因分析 #

报错length(lower) == 1 not TRUE是因为外层integrate调用func2时，会传入向量化的c参数（即c是一个向量），但内层integrate(func1, ...)要求lower必须是标量，导致参数类型冲突。另外代码中存在未定义变量cen.p，实际应为定义好的p。

修正后的数值积分代码 #

alpha = 1
lambda = 4
p = 0.10

# 联合概率密度函数
func1 <- function(t, c, beta) {
  alpha * lambda * exp(-lambda * t^alpha) * beta * exp(-beta * c)
}

# 内层对t积分，用Vectorize支持向量化输入的c
func2 <- Vectorize(function(c, beta) {
  integrate(func1, lower = c, upper = Inf, c = c, beta = beta)$value
})

# 外层对c积分，计算与p的差值
func3 <- function(beta) {
  integrate(func2, lower = 0, upper = Inf, beta = beta)$value - p
}

# 求解beta值
result <- uniroot(func3, lower = 0.001, upper = 10, extendInt = "yes")$root
cat("计算得到的beta值：", round(result, 3), "\n")

运行后输出约0.444，符合预期。

更高效的解析解方式 #

当alpha=1时，Weibull分布退化为指数分布，可手动化简二重积分：

1. 内层对t积分：$\int_{c}^{\infty} \lambda e^{-\lambda t} dt = e^{-\lambda c}$
2. 外层对c积分：$\int_{0}^{\infty} \beta e^{-\beta c} \cdot e^{-\lambda c} dc = \frac{\beta}{\beta + \lambda}$
3. 令积分结果等于p=0.1，解方程$\frac{\beta}{\beta + 4} = 0.1$，得$\beta = \frac{4 \times 0.1}{1 - 0.1} \approx 0.444$

对应的R代码：

lambda = 4
p = 0.10

# 定义解析方程
func_analytic <- function(beta) {
  beta/(beta + lambda) - p
}

# 求解
result_analytic <- uniroot(func_analytic, lower = 0.001, upper = 10)$root
cat("解析解得到的beta值：", round(result_analytic, 3), "\n")

内容的提问来源于stack exchange，提问作者Parviz Shahmirzalou