截断新广义Poisson Lindley分布参数估计：β1估计值异常

针对你在截断新广义Poisson Lindley分布参数MLE估计中beta1结果异常的问题，以下是代码中的关键问题排查及修正：

1. 似然函数未使用传入的数据集参数 #

你的lik函数定义了data参数，但内部直接调用全局变量Tx，逻辑不严谨且易引发错误。修正后让函数依赖传入参数而非全局变量：

lik=function(par, data, Tr){
  alpha1=par[1]; beta1=par[2]
  ll=-sum(log(UT(data, alpha1, beta1, Tr)))
  return(ll)
}

2. 未施加参数非负约束 #

alpha1和beta1均为分布的正参数，默认的Nelder-Mead算法不支持约束，易出现参数非正的迭代结果，导致对数似然计算异常。改用支持约束的L-BFGS-B算法，并设置参数下界：

pp=optim(par=c(alpha1=0.1, beta1=2), 
         fn=lik, 
         data=Tx, 
         Tr=Tr,
         method="L-BFGS-B",
         lower=c(1e-4, 1e-4),  # 避免参数为0或负数
         control = list(parscale=c(alpha1=0.1, beta1=1)))

同时将初始值设为更接近真实值（alpha1=0.1），帮助算法收敛到正确解。

3. 数据生成时的概率归一化 #

原始数据生成中，marginal返回的概率可能因数值精度问题总和略偏离1，需归一化后再采样：

for(i in 1:10000){
  a=0:1000
  pr1=marginal(alpha1,beta1,a)
  pr1=pr1/sum(pr1)  # 归一化概率向量
  y1=sample(a,1,replace=TRUE,prob=pr1)
  t1[i]=y1
}

4. 截断分布的生存函数验证 #

检查UT函数中cf（即CDF(Tr)）的公式是否准确，可通过累加PMF验证：

# 对比公式计算与累加PMF得到的CDF
cdf_check = sum(marginal(alpha1, beta1, 0:Tr))
cf_original = ((alpha1+beta1)*(1+alpha1)^(Tr+2)-(alpha1+alpha1^2+2*alpha1*beta1+alpha1*beta1*Tr+beta1))/((alpha1+beta1)*(1+alpha1)^(Tr+2))
cat("累加PMF的CDF:", cdf_check, "\n公式计算的CDF:", cf_original, "\n")

若两者差异较大，说明cf公式推导有误，需重新核对截断分布的推导过程。

修正后关键代码片段 #

# 修正后的似然函数
lik=function(par, data, Tr){
  alpha1=par[1]; beta1=par[2]
  if(alpha1 <= 0 || beta1 <=0) return(1e10)  # 无效参数返回极大值
  ll=-sum(log(UT(data, alpha1, beta1, Tr)))
  return(ll)
}

# 带约束的MLE估计
pp=optim(par=c(alpha1=0.1, beta1=2), 
         fn=lik, 
         data=Tx, 
         Tr=Tr,
         method="L-BFGS-B",
         lower=c(1e-4, 1e-4),
         control = list(parscale=c(alpha1=0.1, beta1=1)))

# 查看估计结果
pp$par

内容的提问来源于stack exchange，提问作者Sanaullah Khan