Я написал некоторый код на R для выполнения вышеупомянутой задачи в заголовке.Мне было интересно, если кто-нибудь мог бы взглянуть на это и посмотреть, есть ли лучший способ сделать это, так как это заняло у меня много времени, чтобы завершить.
Я опубликую несколько ссылок на несколько страниц из учебника, который я использую для этого, в учебнике уже есть вектор оценок и матрица Гессе, поэтому нет необходимости делать исчисление:
Алгоритм Ньютауна-Рафсона (A.6)
Вектор оценки и матрица Гессиана
Попытка получить этот вывод (2-я таблица)
Вот код, который я написал:
# 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)
Мой вывод
Мое значение для H_alphalambda (3-й ряд, 8-й/ последний столбец) отличается от описанного в учебнике (3-я ссылка), скорее всего, это ошибка автора учебника.Я всегда ищу способы улучшить свои навыки кодирования / эффективность, любые предложения о лучшем способе получить результат были бы очень благодарны.Заранее спасибо.