使用constrOptim实现参数和为1的约束优化时初始值合规却报错求助

问题原因 #

constrOptim的约束逻辑是ui %*% theta - ci >= 0，你只设置了ui=rbind(rep(1,3)), ci=c(1)，这代表的约束是参数之和 >= 1，但你需要的是参数之和等于1的等式约束。虽然初始值0.5+0.3+0.2=1恰好满足sum=1，但函数内部的约束检查可能因浮点精度问题，或者默认的约束容忍度，将其判定为不满足；更关键的是，单一的>=约束无法实现等式约束，必须同时加上<=的约束才能锁定sum=1。

修复代码 #

将约束矩阵ui设置为两行，分别对应sum(theta)>=1和sum(theta)<=1（后者等价于sum(-theta)>= -1），以此实现等式约束：

# Generate example data
set.seed(123)
x1 <- rnorm(1000)
x2 <- rnorm(1000)
y <- rbinom(1000, 1, plogis(0.5 + 0.3*x1 + 0.2*x2)) 

# log-likelihood 
logl <- function(theta,x,y){
  x <- as.matrix(x)
  beta <- theta[1:ncol(x)]
  loglik <- sum(-y*log(1 + exp(-(x%*%beta))) - (1-y)*log(1 + exp(x%*%beta)))
  return(-loglik)
}

# 修改约束：同时设置sum(theta)>=1和sum(theta)<=1
constrOptim(
  theta = c(0.5, 0.3, 0.2), 
  f = logl, 
  x = cbind(1, x1, x2), 
  y = y, 
  hessian = TRUE, 
  method = "L-BFGS-B", 
  control = list(fnscale=1), 
  ui = rbind(rep(1, 3), rep(-1, 3)),  # 两行约束：sum>=1 和 sum<=1
  ci = c(1, -1)                       # 对应约束的阈值
)

额外提示 #

• 原代码中y <- y属于冗余代码，可直接删除，不影响功能运行。
• constrOptim原生仅支持线性不等式约束，若需实现等式约束，必须通过一对方向相反的不等式约束组合来完成。

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