R语言optim()函数初始参数评估失败、参数估计报错排查

通过指定alpha0、alpha1、beta1三个参数，基于自依赖条件泊松分布递推生成模拟数据：设置单条样本长度为200，共生成100组重复样本，初始模拟代码如下：

a0<-5
a1<-0.9
b1<-0.2

l<-rep(1,200)
xs<-rep(0,200)
y<-rep(0,200)
s<-matrix(nrow=100, ncol=200)


xs[1]<-0
l[1]<-1

for (j in 1: 100){
  for (i in 2: 50)
  {
    l[i]<-a0+a1*xs[i-1]+b1*l[i-1]
    xs[i]<-rpois(1,lambda = l[i])
  }
  s[j,1:200]<-xs
}

后续计划调用optim()函数最小化负对数似然，计算100组样本对应的三个参数估计值及对应标准差，但运行失败。初始构造的负对数似然函数与参数估计代码如下：

loglik<-function(theta,x)
{  
  alpha0<-theta[1];
  alpha1<-theta[2];
  beta1<-theta[3]
  
  #lambda
  T<-length(x);
  
  lambda<-rep(1,T);
  likeli<-rep(1,T);
  
  for(t in (2:T))
  {
    
    lambda[t]<-alpha0+alpha1*x[t-1]+beta1*lambda[t-1]; 
    likeli[t]<-((lambda[t]^x[t])*exp(-lambda[t]))/factorial(x[t])
  }
  
  return(-log(prod(likeli)))
}

estimates<-matrix(nrow=100, ncol=3)

for(i in 1:100){
  initial<-c(6,0.8,1)
  res<-optim(initial,loglik,x=s[i,],control=list(maxit=10000),hessian=T)
  estimates[i,] <- res$par
  
}

mean(estimates[,1])
sqrt(var(estimates[,1]))
mean(estimates[,2])
sqrt(var(estimates[,2]))
mean(estimates[,3])
sqrt(var(estimates[,3]))

运行时抛出如下错误，提示无法在初始参数处评估函数：

Error in optim(initial, loglik, x = s[i, ], control = list(maxit = 10000),  : 
  function cannot be evaluated at initial parameters

• 模拟数据生成逻辑存在漏洞：内层生成观测值的循环仅迭代到i=50，单条样本第51位到200位全部保留初始赋值的0，长序列尾部连续0值会导致lambda递推结果异常。
• 似然函数缺少参数合法性校验：泊松分布的强度参数lambda必须严格大于0，初始参数中beta1=1远高于真实值0.2，优化过程中如果探索到导致lambda非正的参数组合，泊松概率计算会返回NaN/非法值，导致函数无法正常求值。
• 似然值计算存在数值下溢：直接对200个0-1区间的概率值做乘积运算，连续相乘会快速下溢为0，后续取对数会得到无意义的-Inf，同样会导致函数评估失败。

1. 修正模拟数据生成逻辑，将内层循环终止值从50改为200，保证200长度的序列全部生成有效观测值。
2. 重写负对数似然函数，提升数值稳定性：
   • 增加参数约束：如果出现alpha0<=0、alpha1<0、beta1<0，或递推过程中任意时刻lambda<=0，直接返回一个极大的惩罚值（如1e10），标记该参数组合不可行。
   • 替换概率乘积取对数的逻辑，直接计算逐点对数泊松概率再求和，从根源避免数值下溢，可直接调用R内置的dpois(..., log=TRUE)计算对数概率。
3. 调整初始参数，将初始beta1值从1调整为更接近真实值的0.3，避免初始点落在非法参数区域。

修正后的完整可运行代码如下：

# 修正后的数据生成代码
a0<-5
a1<-0.9
b1<-0.2

l<-rep(1,200)
xs<-rep(0,200)
s<-matrix(nrow=100, ncol=200)

xs[1]<-0
l[1]<-1

for (j in 1: 100){
  # 循环终止值改为200，生成完整序列
  for (i in 2: 200)
  {
    l[i]<-a0+a1*xs[i-1]+b1*l[i-1]
    xs[i]<-rpois(1,lambda = l[i])
  }
  s[j,1:200]<-xs
}

# 修正后的负对数似然函数
loglik<-function(theta,x)
{  
  alpha0<-theta[1]
  alpha1<-theta[2]
  beta1<-theta[3]
  
  # 参数合法性校验，不满足约束返回大惩罚值
  if(alpha0 <= 0 || alpha1 <0 || beta1 <0){
    return(1e10)
  }
  
  T<-length(x)
  lambda<-rep(1,T)
  ll<-0 # 直接累加对数似然
  
  for(t in 2:T)
  {
    lambda[t]<-alpha0+alpha1*x[t-1]+beta1*lambda[t-1]
    # 递推出非正lambda直接返回惩罚值
    if(lambda[t] <= 0){
      return(1e10)
    }
    # 直接计算对数概率累加，避免下溢
    ll <- ll + dpois(x[t], lambda = lambda[t], log = TRUE)
  }
  
  return(-ll)
}

# 参数估计
estimates<-matrix(nrow=100, ncol=3)
ses<-matrix(nrow=100, ncol=3)

for(i in 1:100){
  # 调整初始参数，beta1设为更接近真实值的0.3
  initial<-c(6,0.8,0.3)
  res<-optim(initial,loglik,x=s[i,],control=list(maxit=10000),hessian=T)
  estimates[i,] <- res$par
  # 可通过海森矩阵逆矩阵的对角线开根计算标准误
  ses[i,] <- sqrt(diag(solve(res$hessian)))
}

# 输出参数估计均值与经验标准差
mean(estimates[,1])
sqrt(var(estimates[,1]))
mean(estimates[,2])
sqrt(var(estimates[,2]))
mean(estimates[,3])
sqrt(var(estimates[,3]))

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