多元Newton-Raphson法参数估计R代码不收敛问题求助

我用R编写了基于多元Newton-Raphson方法估计alpha1、alpha2、beta1、beta2参数的代码，但迭代至1000次仍未收敛，且误差持续增大。已检查数据和函数逻辑无问题，请求排查代码问题，或提供其他参数估计方法。代码如下：

library(base)
library(pracma)
library(readxl)

x <- asuransiOhlssonaft$freq.claim
z <- asuransiOhlssonaft$z
t<-length(x)
t_seq <- seq(1, t, 1)
x_bar <- mean(x)

# first derivative
grad_f <- function(x_0, x, x_bar, t, t_seq, z) {
  alpha_1 <- x_0[1]
  alpha_2 <- x_0[2]
  beta_2 <- x_0[3]
  d_alpha_1 <- t * (log(alpha_1/x_bar) + 1 - digamma(alpha_1) - log((alpha_1/x_bar)+1) - 1) + sum((digamma(alpha_1 + x[t_seq])))
  d_alpha_2 <- t * (-digamma(alpha_2) + digamma(alpha_2 + beta_2)) + sum((digamma(z[t_seq] + alpha_2))-digamma(alpha_2 + beta_2 + x[t_seq]))
  d_beta_2 <- t * (-digamma(beta_2) + digamma(alpha_2 + beta_2)) + sum(digamma(x[t_seq] + beta_2 - z[t_seq])-digamma(alpha_2 + beta_2 + x[t_seq]))
  m <- matrix(c(d_alpha_1, d_alpha_2, d_beta_2),nrow=3,ncol = 1)
  # print(m)
  return (m)
}

# second derivative
hess_f<- function(x_0, x, x_bar, t, t_seq, z) {
  alpha_1 <- x_0[1]
  alpha_2 <- x_0[2]
  beta_2 <- x_0[3]
  f1a1<-t * ((1/alpha_1)-trigamma(alpha_1)-(1/alpha_1)) + sum(trigamma(alpha_1 + x[t_seq]))
  f1a2<-0
  f1b2<-0
  f2a1<-0
  f2a2<-t * (-trigamma(alpha_2) + trigamma(alpha_2 + beta_2)) + sum(trigamma(z[t_seq] + alpha_2)-trigamma(alpha_2 + beta_2 + x[t_seq]))
  f2b2<-t * trigamma(alpha_2+beta_2) - sum(trigamma(alpha_2+beta_2+x[t_seq]))
  f3a1<-0
  f3a2<-t * trigamma(alpha_2+beta_2) - sum(trigamma(alpha_2+beta_2+x[t_seq]))
  f3b2<-t * (-trigamma(beta_2) + trigamma(alpha_2 + beta_2)) + sum(trigamma(x[t_seq] + beta_2 - z[t_seq])-trigamma(alpha_2 + beta_2 + x[t_seq]))
  matrix(c(f1a1,f1a2,f1b2,f2a1,f2a2,f2b2,f3a1,f3a2,f3b2),ncol=3,byrow=TRUE)
}

x <- asuransiOhlssonaft$freq.claim
z <- asuransiOhlssonaft$z
t<-length(x)
t_sequence <- seq(1, t, 1)
x_bar <- mean(x)

x0 <- c(1, 2, 3)
x_old<-x0

# Implement the Newton-Raphson method
tol <- 0.1
max_iter <- 20
iter <- 0
x_lama <- matrix(x0)
while (iter < max_iter) {
  iter <- iter + 1
  
  print(iter)
  print("\n")
  
  x_new <- x_lama - (solve(hess_f(x0, x, x_bar, t, t_sequence, z)) %*% grad_f(x0, x, x_bar, t, t_sequence, z))
  
  print("x new")
  print(abs(x_new - x_old))
  
  if (max(abs(x_new - x_old)) < tol) {
    break
  }
  x_lama <- x_new
  beta_1 <- x_new[1] / x_bar
}

# Print the results
if (iter == max_iter) {
  cat("The Newton-Raphson method did not converge.\n")
} else {
  return(list(x_new[1], x_new[2], beta_1, x_new[3]))
}

代码问题排查 #

• 迭代核心变量错误：循环中始终用初始值x0调用hess_f和grad_f，而非每次迭代更新后的x_lama（当前参数值），导致优化方向完全偏离，参数无法收敛甚至发散。
• 收敛判断失效：x_old在循环中从未更新，导致收敛判断的差值计算完全错误，无法触发终止条件。
• 收敛阈值设置过松：tol=0.1的阈值太大，无法准确判断参数是否收敛。
• 冗余变量定义：重复定义x、z等变量，降低代码可读性。

修复后的Newton-Raphson代码 #

library(pracma)
library(readxl)

# 补充数据加载步骤
asuransiOhlssonaft <- read_excel("你的数据文件路径.xlsx")
x <- asuransiOhlssonaft$freq.claim
z <- asuransiOhlssonaft$z
t <- length(x)
t_sequence <- seq(1, t, 1)
x_bar <- mean(x)

# 一阶导数（梯度）
grad_f <- function(x_0, x, x_bar, t, t_seq, z) {
  alpha_1 <- x_0[1]
  alpha_2 <- x_0[2]
  beta_2 <- x_0[3]
  d_alpha_1 <- t * (log(alpha_1/x_bar) + 1 - digamma(alpha_1) - log((alpha_1/x_bar)+1) - 1) + sum(digamma(alpha_1 + x[t_seq]))
  d_alpha_2 <- t * (-digamma(alpha_2) + digamma(alpha_2 + beta_2)) + sum(digamma(z[t_seq] + alpha_2) - digamma(alpha_2 + beta_2 + x[t_seq]))
  d_beta_2 <- t * (-digamma(beta_2) + digamma(alpha_2 + beta_2)) + sum(digamma(x[t_seq] + beta_2 - z[t_seq]) - digamma(alpha_2 + beta_2 + x[t_seq]))
  return(matrix(c(d_alpha_1, d_alpha_2, d_beta_2), nrow=3, ncol=1))
}

# 二阶导数（Hessian矩阵）
hess_f <- function(x_0, x, x_bar, t, t_seq, z) {
  alpha_1 <- x_0[1]
  alpha_2 <- x_0[2]
  beta_2 <- x_0[3]
  f1a1 <- t * ((1/alpha_1) - trigamma(alpha_1) - (1/alpha_1)) + sum(trigamma(alpha_1 + x[t_seq]))
  f1a2 <- 0
  f1b2 <- 0
  f2a1 <- 0
  f2a2 <- t * (-trigamma(alpha_2) + trigamma(alpha_2 + beta_2)) + sum(trigamma(z[t_seq] + alpha_2) - trigamma(alpha_2 + beta_2 + x[t_seq]))
  f2b2 <- t * trigamma(alpha_2+beta_2) - sum(trigamma(alpha_2+beta_2+x[t_seq]))
  f3a1 <- 0
  f3a2 <- t * trigamma(alpha_2+beta_2) - sum(trigamma(alpha_2+beta_2+x[t_seq]))
  f3b2 <- t * (-trigamma(beta_2) + trigamma(alpha_2 + beta_2)) + sum(trigamma(x[t_seq] + beta_2 - z[t_seq]) - trigamma(alpha_2 + beta_2 + x[t_seq]))
  return(matrix(c(f1a1,f1a2,f1b2,f2a1,f2a2,f2b2,f3a1,f3a2,f3b2), ncol=3, byrow=TRUE))
}

# 初始参数设置
x0 <- c(1, 2, 3)
x_lama <- matrix(x0)
tol <- 1e-6  # 更严格的收敛阈值
max_iter <- 100
iter <- 0
converged <- FALSE

# Newton-Raphson迭代逻辑
while (iter < max_iter) {
  iter <- iter + 1
  # 使用当前参数计算梯度和Hessian
  current_grad <- grad_f(x_lama, x, x_bar, t, t_sequence, z)
  current_hess <- hess_f(x_lama, x, x_bar, t, t_sequence, z)
  
  # 添加小扰动避免Hessian矩阵奇异
  current_hess <- current_hess + diag(3)*1e-8
  x_new <- x_lama - solve(current_hess) %*% current_grad
  
  # 检查收敛条件
  if (max(abs(x_new - x_lama)) < tol) {
    converged <- TRUE
    break
  }
  x_lama <- x_new
}

# 计算beta1
beta_1 <- x_lama[1] / x_bar

# 输出结果
if (converged) {
  cat("收敛成功！\n")
  cat(sprintf("alpha1: %.4f\n", x_lama[1]))
  cat(sprintf("alpha2: %.4f\n", x_lama[2]))
  cat(sprintf("beta1: %.4f\n", beta_1))
  cat(sprintf("beta2: %.4f\n", x_lama[3]))
} else {
  cat("迭代未收敛\n")
}

替代参数估计方法 #

如果修复后的Newton-Raphson仍不收敛，可尝试以下更稳定的方法：

• 拟牛顿法（BFGS）：无需手动计算Hessian矩阵，用近似矩阵替代，稳定性更强，可通过optim函数实现：

  # 定义负对数似然函数（需匹配你的模型逻辑）
  neg_log_lik <- function(params, x, z, x_bar, t) {
    alpha1 <- params[1]
    alpha2 <- params[2]
    beta2 <- params[3]
    # 替换为你的模型负对数似然计算
    ll <- t * (lgamma(alpha1) - alpha1*log(alpha1/x_bar + 1) + alpha1) + sum(lgamma(alpha1 + x) - lgamma(alpha1)) +
          t * (lgamma(alpha2 + beta2) - lgamma(alpha2) - lgamma(beta2)) + sum(lgamma(z + alpha2) + lgamma(x - z + beta2) - lgamma(x + alpha2 + beta2))
    return(-ll)
  }
  
  # 使用BFGS优化
  result <- optim(par = c(1,2,3), fn = neg_log_lik, method = "BFGS", 
                  x = x, z = z, x_bar = x_bar, t = t)
  beta1 <- result$par[1]/x_bar
  cat("alpha1:", result$par[1], "\nalpha2:", result$par[2], "\nbeta1:", beta1, "\nbeta2:", result$par[3], "\n")
  

• EM算法：针对计数数据模型（如索赔频率模型），EM算法通常比Newton-Raphson更稳定，适合含隐变量的场景。
• Nelder-Mead单纯形法：无导数优化方法，适合梯度计算复杂或数值不稳定的情况，同样通过optim函数调用method = "Nelder-Mead"即可。

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