使用R的nlm函数拟合Gompertz模型时beta参数估计异常问题

你的代码存在几个关键问题，导致nlm返回不合理的β值（-999）和极大MSE：

1. 未对参数β施加合理约束 #

Gompertz模型描述体重增长时，β必须为正数：当β>0时，随着时间t增大，exp(-β*t)趋近于0，体重会从初始值p0=40逐渐趋近于上限α=2500，符合鸡的体重增长规律。若β为负，exp(-β*t)会指数级增大，结合log(p0/α)的负值（40<2500），预测体重会趋近于0，完全不符合实际，此时MSE自然极大。你的代码没有限制β的取值范围，nlm可能收敛到这个无意义的区域。

2. 目标函数依赖全局变量（潜在风险） #

mean_squared_error直接调用全局的chicken$weight，虽然当前环境下能运行，但这种写法不健壮，且不利于调试和复用，建议将真实体重作为参数传入函数。

3. 冗余代码未利用 #

你定义了params <- seq(from = 0, to = 1, by = 0.1)但未使用，这部分可以删除。

修正后的代码 #

1. 保留原Gompertz模型（公式正确） #

gompertz_model <- function(times, beta) {
  p0 <- 40
  alpha <- 2500
  
  alpha * exp(log(p0 / alpha) * exp(-beta * times))
}

2. 重构MSE函数（避免全局变量） #

mean_squared_error <- function(beta, model, times, true_weights) {
  predicted_weights <- model(times, beta)
  mean((true_weights - predicted_weights)^2)
}

3. 带约束的参数优化函数 #

给β添加非负约束，当β≤0时返回极大值，引导nlm在合理区域搜索：

find_best_param <- function(model, initial_beta, times, true_weights) {
  # 约束β必须为正，否则返回极大MSE
  constrained_mse <- function(beta) {
    if (beta <= 0) return(1e10)
    mean_squared_error(beta, model, times, true_weights)
  }
  
  res <- nlm(constrained_mse, p = initial_beta)
  return(list(optimal_beta = res$estimate, min_mse = res$minimum))
}

4. 调用优化函数 #

# 假设chicken数据框包含time和weight列
optim_result <- find_best_param(gompertz_model, 
                                initial_beta = 0.1,  # 选择合理的初始值，可尝试0.05-0.2区间
                                times = chicken$time, 
                                true_weights = chicken$weight)

# 查看结果
optim_result$optimal_beta
optim_result$min_mse

额外提示 #

非线性优化对初始值敏感，若一次优化结果不理想，可尝试多个初始值（比如0.05、0.1、0.2），选择MSE最小的结果。

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