如何在mgcv中正确为单调递增GAM模型指定0.1-1.0的上下界？

原mgcv::pcls()文档中的单调递增GAM示例代码如下：

set.seed(1234) # 保证结果可复现
x <- runif(100)*4-1;x <- sort(x);
f <- exp(4*x)/(1+exp(4*x)); y <- f+rnorm(100)*0.1; plot(x,y)
dat <- data.frame(x=x,y=y)
# 无约束样条拟合
f.ug <- gam(y~s(x,k=10,bs="cr")); lines(x,fitted(f.ug))
# 构建单调样条的设计矩阵与约束
sm <- smoothCon(s(x,k=10,bs="cr"),dat,knots=NULL)[[1]]
F <- mono.con(sm$xp);   # 获取单调约束
G <- list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,p=sm$xp,y=y,w=y*0+1)
G$Ain <- F$A;G$bin <- F$b;G$S <- sm$S;G$off <- 0

p <- pcls(G);  # 拟合单调样条

fv<-Predict.matrix(sm,data.frame(x=x))%*%p
lines(x,fv,col=2)

问题分析 #

直接在mono.con()中设置lower=0.1和upper=1.0会报错initial parameters not feasible，原因是：

• mono.con的lower/upper是对样条系数的约束，而非拟合值的约束
• 无约束拟合得到的初始参数sm$xp不满足新增的上下界约束，导致优化无法启动

解决方案：针对拟合值添加上下界约束 #

我们需要直接对拟合值X %*% p构建约束，同时合并单调约束，再调整初始参数到可行域：

修改后的完整代码

set.seed(1234)
x <- runif(100)*4-1;x <- sort(x);
f <- exp(4*x)/(1+exp(4*x)); y <- f+rnorm(100)*0.1; plot(x,y)
dat <- data.frame(x=x,y=y)

# 1. 先做无约束拟合，获取初始平滑参数
f.ug <- gam(y~s(x,k=10,bs="cr")); lines(x,fitted(f.ug), col = "gray")
sm <- smoothCon(s(x,k=10,bs="cr"),dat,knots=NULL)[[1]]

# 2. 构建三种约束：单调递增、拟合值下界0.1、上界1.0
# 单调约束
mono_const <- mono.con(sm$xp)
# 拟合值下界约束：X %*% p >= 0.1 → 等价于 X %*% p - 0.1 >=0
lower_const <- list(A = sm$X, b = rep(0.1, nrow(sm$X)))
# 拟合值上界约束：X %*% p <=1.0 → 等价于 -X %*% p >= -1.0
upper_const <- list(A = -sm$X, b = rep(-1.0, nrow(sm$X)))

# 合并所有约束矩阵和向量
G_Ain <- rbind(mono_const$A, lower_const$A, upper_const$A)
G_bin <- c(mono_const$b, lower_const$b, upper_const$b)

# 3. 调整初始参数到可行域：先把无约束拟合值截断到0.1-1.0，再反推初始p
init_fv <- fitted(f.ug)
init_fv_clamped <- pmax(pmin(init_fv, 1.0), 0.1)
# 用最小二乘反推可行的初始参数p
init_p <- solve(t(sm$X) %*% sm$X, t(sm$X) %*% init_fv_clamped)

# 4. 构建pcls所需的G列表
G <- list(
  X = sm$X,
  C = matrix(0,0,0),  # 无等式约束
  sp = f.ug$sp,       # 用无约束拟合的平滑参数
  p = init_p,         # 用调整后的可行初始参数
  y = y,
  w = rep(1, length(y)),
  Ain = G_Ain,
  bin = G_bin,
  S = sm$S,
  off = 0
)

# 5. 拟合约束模型
p <- pcls(G)

# 6. 计算拟合值并绘图
fv <- Predict.matrix(sm, data.frame(x=x)) %*% p
lines(x, fv, col=2, lwd=2)
# 添加上下界参考线
abline(h=0.1, col="blue", lty=2)
abline(h=1.0, col="blue", lty=2)

关键说明

• 直接针对拟合值X%*%p构建约束，而非对样条系数，这样能精准控制拟合结果的范围
• 初始参数必须调整到满足所有约束的可行域内，这里通过截断无约束拟合值再反推参数实现
• 合并单调约束与上下界约束时，要注意约束的方向（上界约束需要取负号转化为>=形式）

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