如何在R中实现带非负参数约束的迭代加权最小二乘法以复现论文结果？

针对你复现Brookmeyer和Liao 1990年论文的需求，以下是可实现迭代加权最小二乘法（带约束优化）的R包及具体适配方案：

可用R包 #

• nnls：专注非负约束最小二乘，轻量高效，适合基础非负约束场景，可手动嵌套迭代加权逻辑
• lsei：支持线性等式/不等式约束的最小二乘求解，能灵活定义论文中可能涉及的复杂约束（如系数和为1、区间约束等）
• optimx：通用优化框架，可自定义加权最小二乘目标函数，搭配算法实现各类约束（如边界约束、非线性约束）

迭代加权流程适配示例 #

你之前尝试的加权非负最小二乘无法适配迭代流程，核心是需要循环更新权重并重新求解约束最小二乘，直到收敛。以下是具体实现代码：

1. 基于nnls的迭代加权非负最小二乘 #

library(nnls)

iterative_weighted_nnls <- function(X, y, max_iter = 100, tol = 1e-6) {
  # 初始化权重与系数
  weights <- rep(1, length(y))
  beta <- nnls(X, y)$x
  prev_beta <- beta + 1  # 触发首次迭代
  
  iter <- 0
  while (iter < max_iter && max(abs(beta - prev_beta)) > tol) {
    prev_beta <- beta
    # 计算残差并更新权重（可替换为论文指定的权重规则）
    residuals <- y - X %*% beta
    weights <- 1 / (abs(residuals) + 1e-8)  # 加小值避免除零错误
    # 加权变换后执行非负最小二乘
    X_weighted <- X * sqrt(weights)
    y_weighted <- y * sqrt(weights)
    beta <- nnls(X_weighted, y_weighted)$x
    iter <- iter + 1
  }
  list(coefficients = beta, iterations = iter, converged = iter < max_iter)
}

2. 基于lsei的复杂约束迭代加权最小二乘 #

如果论文涉及非负以外的约束（如系数和为1），可使用lsei定义约束条件：

library(lsei)

iterative_weighted_constrained_ls <- function(X, y, max_iter = 100, tol = 1e-6) {
  weights <- rep(1, length(y))
  beta <- rep(1/ncol(X), ncol(X))  # 初始系数设为均匀分布
  prev_beta <- beta + 1
  
  iter <- 0
  while (iter < max_iter && max(abs(beta - prev_beta)) > tol) {
    prev_beta <- beta
    residuals <- y - X %*% beta
    weights <- 1 / (abs(residuals) + 1e-8)
    # 加权变换
    X_weighted <- X * sqrt(weights)
    y_weighted <- y * sqrt(weights)
    
    # 定义约束：beta >=0 且 sum(beta) = 1（可根据论文调整）
    A_eq <- matrix(1, nrow = 1, ncol = ncol(X))  # 等式约束矩阵
    b_eq <- 1  # 等式约束目标值
    A_ineq <- diag(ncol(X))  # 不等式约束矩阵（非负）
    b_ineq <- rep(0, ncol(X))  # 不等式约束下限
    
    # 求解带约束的加权最小二乘
    opt_result <- lsei(X_weighted, y_weighted, A_ineq, b_ineq, A_eq, b_eq)
    beta <- opt_result$X
    iter <- iter + 1
  }
  list(coefficients = beta, iterations = iter, converged = iter < max_iter)
}

3. 基于optimx的自定义约束迭代加权 #

若需要更灵活的约束（如非线性约束），可自定义目标函数搭配optimx的约束算法：

library(optimx)

# 定义加权最小二乘目标函数
weighted_ls_obj <- function(beta, X, y, weights) {
  sum(weights * (y - X %*% beta)^2)
}

iterative_weighted_optimx <- function(X, y, max_iter = 100, tol = 1e-6) {
  weights <- rep(1, length(y))
  beta <- rep(0, ncol(X))
  prev_beta <- beta + 1
  
  iter <- 0
  while (iter < max_iter && max(abs(beta - prev_beta)) > tol) {
    prev_beta <- beta
    residuals <- y - X %*% beta
    weights <- 1 / (abs(residuals) + 1e-8)
    
    # 使用L-BFGS-B算法支持非负边界约束
    opt_res <- optimx(par = beta, fn = weighted_ls_obj, 
                      X = X, y = y, weights = weights,
                      method = "L-BFGS-B", lower = rep(0, ncol(X)))
    beta <- as.numeric(opt_res[1, 1:ncol(X)])
    iter <- iter + 1
  }
  list(coefficients = beta, iterations = iter, converged = iter < max_iter)
}

关键注意事项 #

• 权重更新逻辑需严格匹配论文中的设定（上述示例用的是残差绝对值倒数，需替换为论文指定的权重规则）
• 若迭代收敛慢，可调整max_iter（最大迭代次数）和tol（收敛阈值）参数
• 验证结果时，建议对比论文中的示例数据或模拟数据，确保实现逻辑与论文一致

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