Newton-Raphson’s method for Weibull distribution

Let \(X_{i} \sim \mbox{Weibull}(\alpha, \beta)\) then,

\[\begin{equation*} f(x_{i}|\alpha, \beta) = \dfrac{\alpha}{\beta}\left(\dfrac{x_{i}}{\beta}\right)^{\alpha -1}\exp\left\{-\left(\dfrac{x_{i}}{\beta}\right)^{\alpha}\right\}, \hspace{0.5cm} \mbox{para} \hspace{0.5cm} x_{i}\geq 0. \end{equation*}\]

In this continuous distribution we have that \(\alpha>0\) is the shape parameter and \(\beta>0\) is the scale parameter.

Suposse that \(X_{1}, \cdots, X_{n}\) be a random sample with Weibull distribution \((\alpha,1)\)

Obtain the maximum likelihood estimate using the Newton-Raphson’s method for the alpha parameter for the Weibull distribution.

library(tidyverse)
set.seed(25)

dados <- read_table("Dados_Weibull.txt", col_names = FALSE)

ggplot(dados, aes(x = X1)) +
  geom_histogram(aes(y = after_stat(density)), 
                 bins = 30, 
                 color = "black", 
                 fill = "grey"
  ) +
  geom_density(alpha = 0.2, 
               linetype = 2, 
               fill = "light blue"
  ) +
  labs(x = "Data", y = "Density")

newton.raphson.w <- function(x.dados, alpha.0 = 10, precisao = 1e-7, n = 100) {
  dlogLikW <- function(y) {
    (length(x.dados) / y) + sum(log(x.dados)) - sum(log(x.dados) * (x.dados^y))
  }
  
  ddlogLikW <- function(z) {
    -(length(x.dados) / z^2) - sum((log(x.dados)^2) * (x.dados^z))
  }
  
  for (i in 1:n) {
    alpha.1 <- alpha.0 - dlogLikW(alpha.0) / ddlogLikW(alpha.0)
    
    if (abs(alpha.1 - alpha.0) < precisao) {
      res <- list(alpha.estimado = alpha.1, n.iter = i)
      return(res)
    }
    
    alpha.0 <- alpha.1
  }
  
  print("with the number of iterations there was no convergence.")
}

dados %>%
  pull() %>%
  newton.raphson.w()

$alpha.estimado
[1] 4.965997

$n.iter
[1] 5

#simulation
n <- seq(10, 1000, 10)

estimativas <- sapply(n, function(size) {
  dat <- rweibull(size, shape = 5, scale = 1)
  newton_result <- newton.raphson.w(dat, 2)
  newton_result$alpha.estimado
})

dados <- tibble(
  estim = estimativas,
  sim = n
)

ggplot(dados, aes(x = sim, y = estim)) + 
  geom_point(color = "blue") +
  geom_line(y = 5, linetype = "dashed") +
  labs(x = "Sample size", 
       y = expression(paste("Estimate for ", alpha))
  )

Therefore, the estimated alpha for the data provided was \(\hat{\alpha}=4.9659\). Moreover, we can observe that in the simulation part, the larger the sample, the algorithm will have the better performance.