如何用R求解多函数组合积分约束下的逆高斯分布参数？

已知总函数为积分形式：

C(x)=∫₀^∞ R(x-t)*G(t)dt

其中：

• R(x-t)的定义：

R <- function(t, x){(x-t)*0.0005782559 - 0.1619717}

• G(t)为逆高斯分布，定义如下：

G(t) = sqrt(a^3/(4*pi*b^2*t^3)) * exp((-a*(t-a)^2)/(4*b^2*t))

注：a为均值年龄，b为宽度，满足：

a = integrate(t*G(t), 0, Inf)$value
b = sqrt(integrate((t-a)^2*G(t), 0, Inf)$value/2)

给定条件：

• x=2012，C(x)=269.64
• 假设b/a=1（即b=a）

需要求解参数a和b的值。

用户尝试的R代码：

G <- function(t,a,b=a){
   a <- integrate(G(t,a,b)*t,0,Inf)$value
   b <- sqrt(integrate((t-a)^2*G(t,a,b),0,Inf)$value/2)
   sqrt(a^3/(4*pi*b^2*t^3))*exp((-a*(t-a)^2)/(4*b^2*t))
}
C <- function(y,t){(y-t)*0.0005782559-0.1619717}
g <- Vectorize(function(a) integrate(function(t) G(t=25,a,b=a)*C(y=2012,t=25),lower=0,upper=Inf)$value-269.64)
uniroot(g,c(.01,500))$root

运行后报错：

Error: C stack usage  15927808 is too close to the limit

错误原因 #

1. 无限递归调用：G函数内部反复调用自身计算积分，导致无限递归，触发栈溢出。
2. 参数逻辑混乱：逆高斯分布的a本身就是分布均值，无需再通过积分计算；b的定义也混淆了分布参数与统计量的关系。
3. 积分变量错误：代码中错误地将t=25固定代入积分，而积分变量t应从0到∞遍历，不是固定值。

正确推导与解法 #

已知b=a，可先简化积分表达式：

1. 展开R(x-t)：R(x-t) = k*(x-t) - c，其中k=0.0005782559，c=0.1619717。
2. 拆分积分：

   C(x) = k*x*∫₀^∞ G(t)dt - k*∫₀^∞ t*G(t)dt - c*∫₀^∞ G(t)dt
   

   由于G(t)是概率密度函数，∫₀^∞ G(t)dt=1；∫₀^∞ t*G(t)dt=a（分布均值），代入后得到：

   C(x) = k*(x - a) - c
   

3. 代入已知条件直接解方程：

   269.64 = 0.0005782559*(2012 - a) - 0.1619717
   

   解此方程即可得到a，再由b=a得到b的值。

代码实现 #

方法1：直接解方程（高效）

# 已知参数
k <- 0.0005782559
c <- 0.1619717
x <- 2012
C_x <- 269.64

# 求解a
a <- x - (C_x + c)/k
b <- a  # 因b/a=1

cat("a =", round(a, 2), "\nb =", round(b, 2), "\n")

方法2：数值积分验证（用于确认结果）

如果需要通过积分验证，需避免递归，直接定义逆高斯密度函数：

# 定义逆高斯分布密度（b=a）
G <- function(t, a) {
  sqrt(a^3/(4*pi*a^2*t^3)) * exp((-a*(t - a)^2)/(4*a^2*t))
}

# 定义R函数
R <- function(x, t) {
  (x - t)*0.0005782559 - 0.1619717
}

# 目标函数：计算积分与269.64的差值
target <- function(a) {
  integral <- integrate(function(t) R(2012, t)*G(t, a), 0, Inf)$value
  integral - 269.64
}

# 数值求解
result <- uniroot(target, c(-300000, 0))$root
cat("a =", round(result, 2), "\nb =", round(result, 2), "\n")

结果说明 #

直接解方程得到的a值约为-267000，数值较大，需确认问题中R函数的系数是否正确，但推导逻辑是严谨的。

内容的提问来源于stack exchange，提问作者H Y