如何用Newton–Raphson算法拟合逻辑增长模型及R实现问题

表2.3 不同时间点杂拟谷盗种群数量统计 #

表2.3给出了不同时间点杂拟谷盗（Tribolium confusum）种群的数量统计[103]，统计包含所有发育阶段的甲虫，且食物供应受到严格控制。

种群增长的基础模型为逻辑增长模型，表达式如下：
$$\frac{dN}{dt}=r N (1-\frac{N}{K})$$
其中$N$为种群规模，$t$为时间，$r$为增长率参数，$K$为环境的种群承载能力参数。该微分方程的解为：
$$N_{t}=f(t)=\frac{K N_{0}}{N_{0}+(K-N_{0})\exp(-rt)}$$

问题要求 #

(a) 采用Newton–Raphson方法最小化模型预测值与观测值之间的平方误差和，将逻辑增长模型拟合到上述甲虫种群数据。

(b) 在许多种群建模应用中，会采用对数正态假设：$\log N_t$独立且服从均值为$\log f(t)$、方差为$\sigma^2$的正态分布。分别使用Gauss–Newton和Newton–Raphson方法求该假设下的极大似然估计（MLEs），给出参数估计的标准误差及参数间的相关系数估计，并进行说明。

遇到的问题 #

尝试使用R语言中的pracma包解决问题(b)但未成功，找不到能通过Newton–Raphson算法拟合逻辑增长模型的R包，使用pracma::newtonsys函数也失败了。

手动实现解决方案 #

由于现成R包没有针对该模型的直接实现，我们可以手动编写Gauss-Newton和Newton-Raphson算法来完成拟合。

1. 准备数据

# 加载观测数据
t <- c(0, 8, 28, 41, 63, 79, 97, 117, 135, 154)
N <- c(2, 47, 192, 256, 768, 896, 1120, 896, 1185, 1024)
log_N <- log(N)  # 对数转换后的观测值

2. 定义核心函数

# 逻辑增长模型的对数预测值
log_f <- function(theta, t) {
  K <- theta[1]
  r <- theta[2]
  N0 <- theta[3]
  log(K * N0 / (N0 + (K - N0) * exp(-r * t)))
}

# 负对数似然函数（用于MLE最小化）
nll <- function(theta, t, log_N) {
  K <- theta[1]
  r <- theta[2]
  N0 <- theta[3]
  sigma2 <- theta[4]
  mu <- log_f(c(K, r, N0), t)
  sum((log_N - mu)^2)/(2*sigma2) + length(log_N)/2 * log(2*pi*sigma2)
}

3. Gauss-Newton算法实现

# 计算雅可比矩阵（预测值对K、r、N0的偏导）
jacobian <- function(theta, t) {
  K <- theta[1]
  r <- theta[2]
  N0 <- theta[3]
  denom <- N0 + (K - N0)*exp(-r*t)
  dmu_dK <- (N0 * exp(-r*t)) / denom
  dmu_dr <- (K * N0 * t * exp(-r*t)) / denom
  dmu_dN0 <- (K * (K - N0)*exp(-r*t)) / (N0 * denom)
  cbind(dmu_dK, dmu_dr, dmu_dN0)
}

# Gauss-Newton拟合函数
gauss_newton <- function(theta_init, t, log_N, tol=1e-6, max_iter=100) {
  theta <- theta_init
  iter <- 0
  while(iter < max_iter) {
    mu <- log_f(theta[1:3], t)
    J <- jacobian(theta[1:3], t)
    residuals <- log_N - mu
    # 更新K、r、N0
    theta_update <- theta[1:3] + solve(t(J) %*% J) %*% t(J) %*% residuals
    # 更新sigma²
    sigma2_new <- sum(residuals^2)/length(log_N)
    # 收敛判断
    if(max(abs(theta_update - theta[1:3])) < tol) {
      theta <- c(theta_update, sigma2_new)
      break
    }
    theta[1:3] <- theta_update
    theta[4] <- sigma2_new
    iter <- iter + 1
  }
  # 计算标准误差和相关系数
  info_mat <- t(J) %*% J / theta[4]
  se <- sqrt(diag(solve(info_mat)))
  cov_mat <- solve(info_mat) * theta[4]
  corr_mat <- cov_mat / outer(se, se)
  list(theta=theta, se=se, corr=corr_mat, iter=iter)
}

# 初始值设置：N0取t=0观测值，K取最大观测值，r取0.05，sigma²初始为1
theta_init <- c(K=1185, r=0.05, N0=2, sigma2=1)
gn_result <- gauss_newton(theta_init, t, log_N)

4. Newton-Raphson算法实现

# 负对数似然的梯度
gradient_nll <- function(theta, t, log_N) {
  K <- theta[1]
  r <- theta[2]
  N0 <- theta[3]
  sigma2 <- theta[4]
  mu <- log_f(c(K, r, N0), t)
  residuals <- log_N - mu
  denom <- N0 + (K - N0)*exp(-r*t)
  
  dmu_dK <- (N0 * exp(-r*t)) / denom
  dmu_dr <- (K * N0 * t * exp(-r*t)) / denom
  dmu_dN0 <- (K * (K - N0)*exp(-r*t)) / (N0 * denom)
  
  d_nll_dK <- sum(residuals * dmu_dK)/sigma2
  d_nll_dr <- sum(residuals * dmu_dr)/sigma2
  d_nll_dN0 <- sum(residuals * dmu_dN0)/sigma2
  d_nll_dsigma2 <- -length(log_N)/(2*sigma2) + sum(residuals^2)/(2*sigma2^2)
  
  c(d_nll_dK, d_nll_dr, d_nll_dN0, d_nll_dsigma2)
}

# 负对数似然的Hessian矩阵
hessian_nll <- function(theta, t, log_N) {
  K <- theta[1]
  r <- theta[2]
  N0 <- theta[3]
  sigma2 <- theta[4]
  mu <- log_f(c(K, r, N0), t)
  residuals <- log_N - mu
  n <- length(log_N)
  
  denom <- N0 + (K - N0)*exp(-r*t)
  exp_rt <- exp(-r*t)
  
  # 一阶偏导
  dmu_dK <- (N0 * exp_rt) / denom
  dmu_dr <- (K * N0 * t * exp_rt) / denom
  dmu_dN0 <- (K * (K - N0)*exp_rt) / (N0 * denom)
  
  # 二阶偏导
  d2mu_dK2 <- -N0^2 * exp_rt^2 / denom^2
  d2mu_dKr <- N0 * t * exp_rt * (N0 + (K - N0)*exp_rt*(1 + r*t)) / denom^2
  d2mu_dKN0 <- -K * exp_rt * (N0 + (K - N0)*exp_rt) / (N0^2 * denom^2)
  d2mu_dr2 <- -K * N0 * t^2 * exp_rt * (N0 + (K - N0)*exp_rt*(1 + r*t)) / denom^2
  d2mu_drN0 <- K * t * exp_rt * (K*N0 + (K - N0)^2*exp_rt) / (N0^2 * denom^2)
  d2mu_dN02 <- -K^2 * (K - N0)^2 * exp_rt^2 / (N0^3 * denom^2)
  
  # 构造Hessian矩阵
  H11 <- sum(dmu_dK^2)/sigma2 + sum(residuals * d2mu_dK2)/sigma2
  H12 <- sum(dmu_dK*dmu_dr)/sigma2 + sum(residuals * d2mu_dKr)/sigma2
  H13 <- sum(dmu_dK*dmu_dN0)/sigma2 + sum(residuals * d2mu_dKN0)/sigma2
  H14 <- sum(residuals * dmu_dK)/sigma2^2
  
  H22 <- sum(dmu_dr^2)/sigma2 + sum(residuals * d2mu_dr2)/sigma2
  H23 <- sum(dmu_dr*dmu_dN0)/sigma2 + sum(residuals * d2mu_drN0)/sigma2
  H24 <- sum(residuals * dmu_dr)/sigma2^2
  
  H33 <- sum(dmu_dN0^2)/sigma2 + sum(residuals * d2mu_dN02)/sigma2
  H34 <- sum(residuals * dmu_dN0)/sigma2^2
  
  H44 <- n/(2*sigma2^2) - sum(residuals^2)/sigma2^3
  
  hess <- matrix(0, nrow=4, ncol=4)
  hess[1,1] <- H11; hess[1,2] <- H12; hess[1,3] <- H13; hess[1,4] <- H14
  hess[2,1] <- H12; hess[2,2] <- H22; hess[2,3] <- H23; hess[2,4] <- H24
  hess[3,1] <- H13; hess[3,2] <- H23; hess[3,3] <- H33; hess[3,4] <- H34
  hess[4,1] <- H14; hess[4,2] <- H24; hess[4,3] <- H34; hess[4,4] <- H44
  
  hess
}

# Newton-Raphson拟合函数
newton_raphson <- function(theta_init, t, log_N, tol=1e-6, max_iter=100) {
  theta <- theta_init
  iter <- 0
  while(iter < max_iter) {
    grad <- gradient_nll(theta, t, log_N)
    hess <- hessian_nll(theta, t, log_N)
    # 参数更新
    theta_update <- theta - solve(hess) %*% grad
    # 收敛判断
    if(max(abs(theta_update - theta)) < tol) {
      theta <- theta_update
      break
    }
    theta <- theta_update
    iter <- iter + 1
  }
  # 计算标准误差和相关系数
  cov_mat <- solve(hessian_nll(theta, t, log_N))
  se <- sqrt(diag(cov_mat))
  corr_mat <- cov_mat / outer(se, se)
  list(theta=theta, se=se, corr=corr_mat, iter=iter)
}

# 运行Newton-Raphson
nr_result <- newton_raphson(theta_init, t, log_N)

结果说明

运行上述代码后，可得到两种算法的参数估计结果，包括$K$、$r$、$N_0$、$\sigma^2$的估计值，以及对应的标准误差和相关系数矩阵：

• 标准误差：反映参数估计的不确定性，数值越小说明估计精度越高。
• 相关系数矩阵：展示参数间的线性相关程度，若两个参数的相关系数接近±1，说明存在强共线性，可能导致估计不稳定。

注意事项：

1. 初始值选择对收敛至关重要，建议$N_0$取t=0的观测值，$K$取最大观测种群数，$r$取较小正数。
2. Gauss-Newton计算更简单，依赖雅可比矩阵近似Hessian，适合非线性最小二乘场景；Newton-Raphson需精确计算Hessian，推导复杂但收敛速度更快。
3. 若算法不收敛，可调整收敛阈值、最大迭代次数或更换初始值。

内容的提问来源于stack exchange，提问作者Ali Mohamed Hassan