R优化问题：alpha1、alpha2、theta估计准确但beta1、beta2估计偏差

我使用R语言的optim()函数通过极大似然估计（MLE）拟合包含alpha1、beta1、alpha2、beta2、theta五个参数的Lindley（BGPL）分布，但MLE结果存在偏差，尤其是beta1和beta2的估计值始终偏离真实值。已尝试调整不同样本量和初始值，问题仍未解决。

以下是完整R代码：

rm(list = ls())

# 真实参数
alpha1 <- 1.5; beta1 <- 3; alpha2 <- 1.5; beta2 <- 3; theta <- 2

# 定义计算m1和m2的函数
m1 <- function(alpha1, beta1) {
  numerator <- alpha1^2 * (beta1 + alpha1 - exp(-1) + 1)
  denominator <- (alpha1 + beta1) * (alpha1 - exp(-1) + 1)^2
  numerator / denominator
}

m2 <- function(alpha2, beta2) {
  numerator <- alpha2^2 * (beta2 + alpha2 - exp(-1) + 1)
  denominator <- (alpha2 + beta2) * (alpha2 - exp(-1) + 1)^2
  numerator / denominator
}

# X1的边缘PMF
marginal <- function(alpha1, beta1, x1) {
  (alpha1^2 * (1 + alpha1 + beta1 + beta1 * x1)) / ((alpha1 + beta1) * (alpha1 + 1)^(x1 + 2))
}

# 给定X1时X2的条件PMF
conditional <- function(alpha1, beta1, alpha2, beta2, x1, x2, theta) {
  f2 <- (alpha2^2 * (1 + alpha2 + beta2 + beta2 * x2)) / ((alpha2 + beta2) * (alpha2 + 1)^(x2 + 2))
  dep <- 1 + theta * (exp(-x1) - m1(alpha1, beta1)) * (exp(-x2) - m2(alpha2, beta2))
  pmf <- f2 * dep
  return(pmf)
}

# 利用BGPL的边缘和条件PMF生成随机数
t1 <- numeric(1000)
t2 <- numeric(1000)

for (i in 1:1000) {
  a <- 0:1000
  pr1 <- marginal(alpha1, beta1, a)
  pr1 <- pr1 / sum(pr1)
  y1 <- sample(a, 1, replace = TRUE, prob = pr1)
  
  b <- 0:1000
  pr2 <- conditional(alpha1, beta1, alpha2, beta2, y1, b, theta)
  y2 <- sample(b, 1, replace = TRUE, prob = pr2)
  
  t1[i] <- y1
  t2[i] <- y2
}

# X ~ GPL(alpha1, beta1)的边缘PMF
GPLx <- function(x, alpha1, beta1) {
  p1 <- alpha1^2 * (1 + alpha1 + beta1 + beta1 * x)
  p2 <- (alpha1 + beta1) * (alpha1 + 1)^(x + 2)
  p1 / p2
}

# Y ~ GPL(alpha2, beta2)的边缘PMF
GPLy <- function(y, alpha2, beta2) {
  p1 <- alpha2^2 * (1 + alpha2 + beta2 + beta2 * y)
  p2 <- (alpha2 + beta2) * (alpha2 + 1)^(y + 2)
  p1 / p2
}

# BGPL的联合PMF
dBGPL <- function(x, y, alpha1, beta1, alpha2, beta2, theta) {
  m1 <- (alpha1^2 * (beta1 + alpha1 - exp(-1) + 1)) /
    ((alpha1 + beta1) * (alpha1 - exp(-1) + 1)^2)
  m2 <- (alpha2^2 * (beta2 + alpha2 - exp(-1) + 1)) /
    ((alpha2 + beta2) * (alpha2 - exp(-1) + 1)^2)
  
  mx <- GPLx(x, alpha1, beta1)
  my <- GPLy(y, alpha2, beta2)
  
  B1 <- mx * my
  B2 <- 1 + theta * (exp(-x) - m1) * (exp(-y) - m2)
  B1 * B2
}

# 负对数似然函数
lik <- function(pars, data) {
  pa1 <- pars[1]
  ba1 <- pars[2]
  pa2 <- pars[3]
  ba2 <- pars[4]
  tha <- pars[5]
  x <- t1
  y <- t2
  
  -sum(log(dBGPL(x, y, pa1, ba1, pa2, ba2, tha)))
}

# 运行优化器
result <- optim(par = c(1, 1, 1, 1, 1), fn = lik, data = c(t1, t2))
print(result)

排查与优化建议 #

1. 参数约束与优化算法调整

   • optim()默认Nelder-Mead算法对带约束的参数拟合效果有限，建议改用带边界约束的L-BFGS-B算法，同时给beta1、beta2等非负参数设置下限：

     result <- optim(par = c(1.5, 3, 1.5, 3, 2), fn = lik, data = list(x=t1, y=t2),
                    method = "L-BFGS-B", lower = c(1e-6, 1e-6, 1e-6, 1e-6, -Inf))
     

   • 修正数据传递问题：原代码data = c(t1,t2)会合并向量，且似然函数依赖全局变量t1/t2，改为传递列表避免全局变量耦合。

2. 联合PMF合理性验证

   • 检查dBGPL中m1、m2的计算是否严格匹配BGPL分布的理论公式，尤其是exp(-1)的推导逻辑。
   • 确保conditional和dBGPL中1 + theta*(...)部分始终为正，避免出现负概率值导致似然函数异常。

3. 样本生成优化

   • 样本生成时0:1000范围过大，可先计算边缘分布的99.9%分位数，截断到合理取值范围，减少无效计算。
   • 对比生成样本的直方图与理论PMF，确认样本生成过程无逻辑错误。

4. 初始值精细化

   • 不要仅用默认初始值，先通过矩估计得到alpha1、beta1等参数的初始近似值，再传入优化器，提升收敛效率。

内容的提问来源于stack exchange，提问作者Abdul Malik