# Using R to maximize a two parameter Weibull model via multivariate extension of Newton-Raphson method

I am just getting back into using R for the first time in a while, and wrote some code to perform the aforementioned task in the title. I was wondering if anyone could take a look at it and see if there is a better way to do this, as this took me a long time to complete.

I'll post some links to a few pages from the textbook that I'm using for this, the textbook already gives the score vector and Hessian matrix so no need to do any calculus:

Newtown-Raphson Algorithm (A.6)

Score Vector and Hessian Matrix

Trying to get this output (2nd table)

Here is the code that I wrote:

# 10 Uncensored Obs. from Weibull model:
weibull <- c(2.57, 0.58, 0.82, 1.02, 0.78, 0.46, 1.04, 0.43, 0.69, 1.37)
alpha <- 1 # Initial guess for alpha parameter
lambda <- 10 / sum(weibull) # Initial guess for lambda parameter
step <- 0 # Counter

# Score vector:
u_lambda <- (10 / lambda) - sum(weibull^alpha)
u_alpha <- (10 / alpha) + sum(log(weibull)) - lambda*sum(weibull*log(weibull))

# Hessian matrix:
H_lambdalambda <- (-10 / lambda^2)
H_alphaalpha <- (-10 / alpha^2) - lambda*sum(weibull^alpha * log(weibull)^2)
H_alphalambda <- -sum(weibull^alpha * log(weibull))

fun1 <- function(alpha){
a <- weibull^alpha
return (sum(a))
}

fun2 <- function(alpha){
a <- weibull^alpha
b <- log(weibull)
return (sum(a*b))
}

fun3 <- function(alpha){
a <- weibull^alpha
b <- log(weibull)^2
return (sum(a*b))
}

for (i in 1:100){
step[i] <- i-1
u_lambda[i] <- ((10 / lambda[i]) - fun1(alpha[i]))
u_alpha[i] <- ((10 / alpha[i]) + sum(log(weibull)) - lambda[i]*fun2(alpha[i]))
H_lambdalambda[i] <- -(10 / lambda[i]^2)
H_alphaalpha[i] <- (-10/alpha[i]^2 - lambda[i]*fun3(alpha[i]))
H_alphalambda[i] <- (-sum(fun2(alpha[i])))

# Convergence of algorithm is declared under this condition:
if ((abs(u_lambda[i]) < 0.1) && (abs(u_alpha[i]) < 0.1)){
break
}

hessian_matrix <- cbind(c(H_lambdalambda[i], H_alphalambda[i]), c(H_alphalambda[i], H_alphaalpha[i]))
score_vector <- cbind(c(u_lambda[i], u_alpha[i]))
parameters <- c(lambda[i], alpha[i]) - solve(hessian_matrix) %*% score_vector
lambda[i+1] <- parameters[1]
alpha[i+1] <- parameters[2]
}

data.frame(step, lambda, alpha, u_lambda, u_alpha, H_lambdalambda, H_alphaalpha, H_alphalambda)


My Output:

My value for H_alphalambda (3rd row, 8th/last column) is different than that from the textbook (3rd link), most likely it's an error made by the textbook author. I am always looking for ways to improve my coding skills/efficiency, any suggestions on a better way to get to the output would be greatly appreciated. Thanks in advance.

• Just for reference, the C style for loop stuff isn't necessary (don't need to declare i or iterate it). So feel free to rid yourself of step. Feb 28 '19 at 18:46