I want to write an R-code to compute the MLEs of (µ; σ2; λ) with Newton-Raphson algorithm and try initial values λ = 0.5、.... I'm trying to run the following code in R, but I'm getting an error I'm not sure what part of the formula is incorrect. Any help would be greatly appreciated.
obs.loglik <- function(y, theta) {
mu <- theta[1]
sigma2 <- theta[2]
lambda <- theta[3]
n <- length(y)
-0.5*n*log(2*pi) - 0.5*n*log(sigma2) - sum((y^lambda - 1)/lambda - mu)^2/(2*sigma2) +
(lambda - 1)*(sum(log(y)))
}
library(maxLik)
mle <- maxLik(obs.loglik, start=c(mean(y), var(y), 0.5))
# Error in fnOrig(theta, ...) :
# argument "theta" is missing, with no default```
If theta is fixed, we could use it as a global variable without directly providing as input.
obs.loglik <- function(y) {
mu <- theta[1]
sigma2 <- theta[2]
lambda <- theta[3]
n <- length(y)
-0.5*n*log(2*pi) - 0.5*n*log(sigma2) - sum((y^lambda - 1)/lambda - mu)^2/(2*sigma2) +
(lambda - 1)*(sum(log(y)))
}
library(maxLik)
y <- 1:10 # example y
theta <- c(mean(y), var(y), 0.5) # Define theta outside
mle <- maxLik(obs.loglik, start=y)
If theta depends on y, we could do this instead.
obs.loglik <- function(y) {
mu <- mean(y) # depends on input
sigma2 <- var(y) # depends on input
lambda <- theta[3]
n <- length(y)
-0.5*n*log(2*pi) - 0.5*n*log(sigma2) - sum((y^lambda - 1)/lambda - mu)^2/(2*sigma2) +
(lambda - 1)*(sum(log(y)))
}