# Newton Raphson Method for Several Variables in R Studio

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

``````library(base)
library(pracma)

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.

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

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.