如何用R的optim()实现带约束的三次样条参数估计

数据老陈

2026-6-15

带约束的三次函数拟合实现方案（基于R语言）

需求说明

现有数据集Rand_Data，需拟合三次函数：

MyFN = function(a, b, c, d) return(a * Rand_Data$x^3 + b * Rand_Data$x^2 + c * Rand_Data$x + d)

要求用R的优化函数估计参数a、b、c、d，目标是最小化拟合值与Rand_Data$Val的绝对差之和，同时满足两个约束：

第100个位置的拟合值小于2：MyFN[100] < 2
拟合值总和与观测值总和的差值绝对值小于0.001：|sum(MyFN) - sum(Rand_Data$Val)| < 0.001

实现方案

方案1：惩罚项法（通用型，支持非线性约束）

通过在目标函数中加入惩罚项，将约束转化为无约束优化问题，适合各类约束场景。

步骤1：定义带惩罚的目标函数

# 目标函数：L1损失+约束惩罚
objective_fn <- function(par) {
  a <- par[1]
  b <- par[2]
  c <- par[3]
  d <- par[4]
  
  # 计算所有拟合值
  fitted_vals <- a * Rand_Data$x^3 + b * Rand_Data$x^2 + c * Rand_Data$x + d
  
  # 基础损失：绝对差之和
  loss <- sum(abs(fitted_vals - Rand_Data$Val))
  
  # 约束1：第100个拟合值不满足<2时，加入线性惩罚
  if (fitted_vals[100] >= 2) {
    loss <- loss + 1e8 * (fitted_vals[100] - 2)
  }
  
  # 约束2：总和差值不满足<0.001时，加入线性惩罚
  sum_diff <- abs(sum(fitted_vals) - sum(Rand_Data$Val))
  if (sum_diff >= 0.001) {
    loss <- loss + 1e8 * (sum_diff - 0.001)
  }
  
  return(loss)
}

步骤2：设置初始参数

建议用普通最小二乘（OLS）拟合的三次模型参数作为初始值，提升优化效率：

# OLS拟合三次模型获取初始参数
ols_fit <- lm(Val ~ poly(x, 3, raw = TRUE), data = Rand_Data)
# 调整参数顺序为a(x³系数)、b(x²系数)、c(x系数)、d(截距)
init_par <- c(coef(ols_fit)[4], coef(ols_fit)[3], coef(ols_fit)[2], coef(ols_fit)[1])

步骤3：调用optim优化

# 执行优化（选用BFGS方法，适用于光滑目标函数）
optim_result <- optim(par = init_par, fn = objective_fn, method = "BFGS")

# 提取最优参数
best_par <- optim_result$par
names(best_par) <- c("a", "b", "c", "d")
print(best_par)

# 验证约束是否满足
fitted_final <- best_par[1] * Rand_Data$x^3 + best_par[2] * Rand_Data$x^2 + best_par[3] * Rand_Data$x + best_par[4]
cat("第100个拟合值：", fitted_final[100], "\n")
cat("拟合值与观测值总和差值：", abs(sum(fitted_final) - sum(Rand_Data$Val)), "\n")

方案2：constrOptim法（线性约束专用）

本次需求中的两个约束均为线性约束，可直接使用R内置的constrOptim()函数，无需手动加惩罚，约束满足更严格。

步骤1：构造线性约束矩阵

将约束转化为A %*% par <= b的形式：

# 提取数据中的固定值
x100 <- Rand_Data$x[100]
n <- nrow(Rand_Data)
sum_x3 <- sum(Rand_Data$x^3)
sum_x2 <- sum(Rand_Data$x^2)
sum_x <- sum(Rand_Data$x)
sum_val <- sum(Rand_Data$Val)
epsilon <- 1e-8  # 将严格小于转为小于等于的极小偏移量

# 约束矩阵A与约束向量b
A <- rbind(
  c(x100^3, x100^2, x100, 1),          # 约束1：a*x100³ + b*x100² + c*x100 + d ≤ 2-ε
  c(sum_x3, sum_x2, sum_x, n),         # 约束2：sum(fitted) ≤ sum_val + 0.001-ε
  c(-sum_x3, -sum_x2, -sum_x, -n)      # 约束3：sum(fitted) ≥ sum_val - 0.001+ε
)

b <- c(
  2 - epsilon,
  sum_val + 0.001 - epsilon,
  -sum_val + 0.001 - epsilon
)

步骤2：定义无约束目标函数

# 仅计算绝对差之和的目标函数
objective_fn_unconstrained <- function(par) {
  a <- par[1]
  b <- par[2]
  c <- par[3]
  d <- par[4]
  fitted_vals <- a * Rand_Data$x^3 + b * Rand_Data$x^2 + c * Rand_Data$x + d
  sum(abs(fitted_vals - Rand_Data$Val))
}

步骤3：调用constrOptim优化

# 执行约束优化（梯度设为NULL，使用数值梯度）
constr_result <- constrOptim(theta = init_par,
                             f = objective_fn_unconstrained,
                             grad = NULL,
                             ui = A,
                             ci = b)

# 提取最优参数
best_par_constr <- constr_result$par
names(best_par_constr) <- c("a", "b", "c", "d")
print(best_par_constr)

# 验证约束
fitted_constr <- best_par_constr[1] * Rand_Data$x^3 + best_par_constr[2] * Rand_Data$x^2 + best_par_constr[3] * Rand_Data$x + best_par_constr[4]
cat("第100个拟合值：", fitted_constr[100], "\n")
cat("拟合值与观测值总和差值：", abs(sum(fitted_constr) - sum_val), "\n")

注意事项

惩罚项法中的惩罚系数（如1e8）需根据数据规模调整：系数太小约束可能失效，太大可能导致优化不稳定。
constrOptim()仅支持线性约束，若后续需求出现非线性约束，建议优先使用惩罚项法。
若优化结果不理想，可尝试更换optim()的方法（如"L-BFGS-B"）或调整初始参数。

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