Newton Raphson Method for Several Variables in R Studio

[thisisdumpppp]

so I made an R code to estimate parameters (alpha1, alpha2, beta1, and beta2). Here's my code:

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]))
}

the base function:
base function: t * (alpha_1 * log(beta_1) - lgamma(alpha_2) - lgamma(beta_2) + lgamma(alpha_2 + beta_2) - lgamma(alpha_1)- (x_bar+alpha_1)*log(beta_1+1))+sum(lgamma(z[t_seq] + alpha_2) + lgamma(x[t_seq] + beta_2 - z[t_seq]) - lgamma(alpha_2 + beta_2 + x[t_seq])+lgamma(alpha_1+x[t_seq]))

I have no idea if there is something wrong with my code, but until 1000 iterations, the result still "The Newton-Raphson method did not converge". When I check the error by iteration, turns out the error is keep increasing, not decreasing. I'm scared that this long iterations is caused by my mistake on my code. Please help me. Or if there other ways to estimate these parameters than using newton raphson?

Thank you so much!

I already checked for the data and the function, but it is already correct.

[technocrat]

Check the data again with VGAMExtra::newtonRaphson.basic(), which is a mature, well exercised implementation.