### Test the 'finalHessian' argument of optimization routines
library(maxLik)
set.seed( 4 )
# log-likelihood function, gradient, and Hessian for 1-parameter case (exponential distribution)
ll1i <- function(theta) {
if(!all(theta > 0))
return(NA)
log(theta) - theta*t
}
ll1 <- function(theta) sum( log(theta) - theta*t )
gr1i <- function(theta) 1/theta - t
gr1 <- function(theta) sum( 1/theta - t )
hs1 <- function(theta) -100/theta^2
t <- rexp( 100, 2 )
## the same functions for 2-variable case (normal distribution)
ll2 <- function( param ) {
## log likelihood function
mu <- param[ 1 ]
sigma <- param[ 2 ]
if(!(sigma > 0))
return(NA)
# to avoid warnings in the output
N <- length( x )
llValue <- -0.5 * N * log( 2 * pi ) - N * log( sigma ) -
0.5 * sum( ( x - mu )^2 / sigma^2 )
return( llValue )
}
## log likelihood function (individual observations)
ll2i <- function( param ) {
mu <- param[ 1 ]
sigma <- param[ 2 ]
if(!(sigma > 0))
return(NA)
# to avoid warnings in the output
llValues <- -0.5 * log( 2 * pi ) - log( sigma ) -
0.5 * ( x - mu )^2 / sigma^2
return( llValues )
}
gr2 <- function( param ) {
## function to calculate analytical gradients
mu <- param[ 1 ]
sigma <- param[ 2 ]
N <- length( x )
llGrad <- c( sum( ( x - mu ) / sigma^2 ),
- N / sigma + sum( ( x - mu )^2 / sigma^3 ) )
return( llGrad )
}
## function to calculate analytical gradients (individual observations)
gr2i <- function( param ) {
mu <- param[ 1 ]
sigma <- param[ 2 ]
llGrads <- cbind( ( x - mu ) / sigma^2,
- 1 / sigma + ( x - mu )^2 / sigma^3 )
return( llGrads )
}
## function to calculate analytical Hessians
hs2 <- function( param ) {
mu <- param[ 1 ]
sigma <- param[ 2 ]
N <- length( x )
llHess <- matrix( c(
N * ( - 1 / sigma^2 ),
sum( - 2 * ( x - mu ) / sigma^3 ),
sum( - 2 * ( x - mu ) / sigma^3 ),
N / sigma^2 + sum( - 3 * ( x - mu )^2 / sigma^4 ) ),
nrow = 2, ncol = 2 )
return( llHess )
}
x <- rnorm(100, 1, 2)
## NR
# Estimate with only function values (single parameter)
a <- maxLik( ll1i, gr1i, start = 1, method = "NR" )
summary(a )
b <- maxLik( ll1i, gr1i, start = 1, method = "NR", finalHessian="bhhh")
# should issue a warning as BHHH not possible
summary(b )
c <- maxLik( ll1i, gr1i, start = 1, method = "NR", finalHessian=FALSE)
summary(c)
## (vector parameter)
a <- maxLik( ll2, gr2, start = c(0,1), method = "NR" )
summary(a )
b <- maxLik( ll2, gr2, start = c(0,1), method = "NR", finalHessian="bhhh")
# should issue a warning as BHHH not possible
summary(b )
c <- maxLik( ll2, gr2, start = c(0,1), method = "NR", finalHessian=FALSE)
summary(c)
## BFGSR
# Estimate with only function values (single parameter)
a <- maxLik( ll1i, gr1i, start = 1, method = "BFGSR" )
summary(a )
b <- maxLik( ll1i, gr1i, start = 1, method = "BFGSR", finalHessian="bhhh")
# should issue a warning as BHHH not possible
summary(b )
c <- maxLik( ll1i, gr1i, start = 1, method = "BFGSR", finalHessian=FALSE)
summary(c)
# Estimate with only function values (vector parameter)
a <- maxLik( ll2, gr2, start = c(0,1), method = "BFGSR" )
summary(a )
b <- maxLik( ll2, gr2, start = c(0,1), method = "BFGSR", finalHessian="bhhh")
# should issue a warning as BHHH not possible
summary(b )
c <- maxLik( ll2, gr2, start = c(0,1), method = "BFGSR", finalHessian=FALSE)
summary(c)
### Nelder-Mead
## Individual observations only
b <- maxLik( ll2i, start = c(0,1), method = "NM", finalHessian="bhhh")
summary(b)
## Individual observations, summed gradient
b <- maxLik( ll2i, gr2, start = c(0,1), method = "NM", finalHessian="bhhh")
# should issue a warning as BHHH not selected
# (yes, could do it based on individual likelihood and numeric gradient)
summary(b)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.