含floor函数的非线性回归遇奇异梯度矩阵错误的解决方法

错误根源 #

你的模型包含floor()阶梯函数，这是不连续且不可导的。而nls（包括minpack.lm中的实现）依赖参数的梯度（导数）进行迭代优化，阶梯函数在断点处梯度为0或不存在，直接导致梯度矩阵奇异，无法完成参数估计。

替代解决方案 #

1. 用连续函数近似阶梯模型 #

把floor()替换成可导的连续近似函数（比如用sigmoid模拟阶梯跳变），让nls能正常计算梯度：

# 定义带连续近似的模型函数
f_approx <- function(x, k, s, t) {
  z <- (x/t)^s
  # 用sigmoid函数模拟阶梯效果，10是控制跳变陡峭程度的参数
  floor_approx <- z - 0.5 + 0.5*tanh(10*(z - round(z)))
  10/(1 + k*floor_approx)
}

# 拟合模型（初始值尽量接近真实生成参数）
start_values <- c(k=1.5, s=0.7, t=2.5)
fit_approx <- nls(y ~ f_approx(x, k, s, t),
                  data = dat,
                  control = nls.control(maxiter = 1000))
summary(fit_approx)

2. 网格搜索+最小二乘法 #

因为只有3个参数，直接遍历参数取值范围，计算残差平方和找最优组合：

# 设置参数搜索范围
k_vals <- seq(1.0, 2.0, by=0.05)
s_vals <- seq(0.5, 0.9, by=0.02)
t_vals <- seq(2.0, 3.0, by=0.05)

# 生成所有参数组合并计算残差平方和
param_grid <- expand.grid(k=k_vals, s=s_vals, t=t_vals)
param_grid$rss <- apply(param_grid, 1, function(params) {
  y_pred <- 10/(1 + params[1]*floor((dat$x/params[3])^params[2]))
  sum((dat$y - y_pred)^2)
})

# 提取最优参数组合
best_params <- param_grid[which.min(param_grid$rss), ]
print(best_params)

3. 非梯度优化方法（optim函数） #

用不需要计算梯度的优化算法（比如Nelder-Mead）直接最小化残差平方和：

# 定义残差平方和计算函数
rss_func <- function(params) {
  k <- params[1]
  s <- params[2]
  t <- params[3]
  y_pred <- 10/(1 + k*floor((dat$x/t)^s))
  sum((dat$y - y_pred)^2)
}

# 初始值贴近真实参数，提升优化效率
start_vals <- c(k=1.5, s=0.7, t=2.5)
fit_optim <- optim(start_vals, rss_func, method="Nelder-Mead")

# 输出最优参数
print(fit_optim$par)

4. 分段拟合 #

由于数据是阶梯状，每个阶梯内y值恒定，可以先找出x的断点，再针对每个阶梯段的x范围反向推导参数：

# 找出y值突变的x断点
break_indices <- which(diff(dat$y) != 0)
x_breaks <- dat$x[break_indices]

# 每个阶梯对应floor((x/t)^s)的一个整数值，可通过每个段的x范围计算参数关系
# 例如，第一个阶梯y=10，对应floor((x/t)^s)=0 → x < t
# 第二个阶梯y=10/(1+1.5*1)≈4，对应floor((x/t)^s)=1 → t ≤ x < t*(2)^(1/s)
# 以此类推，可通过断点x值建立方程求解k,s,t

重要提示 #

不管用哪种方法，初始参数尽量贴近真实生成值（比如你用的k=1.5、s=0.7、t=2.5），能大幅提升拟合效率和准确性。

内容的提问来源于stack exchange，提问作者Juan Pablo Molano Gallardo