Hessian matrix

This function extracts the Hessian of the objective function at optimum. The Hessian information should be supplied by the underlying optimization algorithm, possibly by an approximation.

hessian(x, ...)
# Default S3 method
hessian(x, ...)

Arguments

x

an optimization result of class ‘maxim’ or ‘maxLik’

...

other arguments for methods

Value

A numeric matrix, the Hessian of the model at the estimated parameter values. If the maximum is flat, the Hessian is singular. In that case you may want to invert only the non-singular part of the matrix. You may also want to fix certain parameters (see activePar).

Author

Ott Toomet

Examples

# log-likelihood for normal density
# a[1] - mean
# a[2] - standard deviation
ll <- function(a) sum(-log(a[2]) - (x - a[1])^2/(2*a[2]^2))
x <- rnorm(100) # sample from standard normal
ml <- maxLik(ll, start=c(1,1))
# ignore eventual warnings "NaNs produced in: log(x)"
summary(ml) # result should be close to c(0,1)
#> --------------------------------------------
#> Maximum Likelihood estimation
#> Newton-Raphson maximisation, 7 iterations
#> Return code 1: gradient close to zero (gradtol)
#> Log-Likelihood: -59.61072 
#> 2  free parameters
#> Estimates:
#>      Estimate Std. error t value Pr(> t)    
#> [1,] -0.03216    0.11009  -0.292    0.77    
#> [2,]  1.10088    0.07784  14.142  <2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> --------------------------------------------
hessian(ml) # How the Hessian looks like
#>           [,1]      [,2]
#> [1,] -82.51533    0.0000
#> [2,]   0.00000 -165.0307
sqrt(-solve(hessian(ml))) # Note: standard deviations are on the diagonal
#>           [,1]       [,2]
#> [1,] 0.1100862 0.00000000
#> [2,] 0.0000000 0.07784266
#
# Now run the same example while fixing a[2] = 1
mlf <- maxLik(ll, start=c(1,1), activePar=c(TRUE, FALSE))
summary(mlf) # first parameter close to 0, the second exactly 1.0
#> --------------------------------------------
#> Maximum Likelihood estimation
#> Newton-Raphson maximisation, 2 iterations
#> Return code 1: gradient close to zero (gradtol)
#> Log-Likelihood: -60.59652 
#> 1  free parameters
#> Estimates:
#>      Estimate Std. error t value Pr(> t)
#> [1,] -0.03216    0.10000  -0.322   0.748
#> [2,]  1.00000    0.00000      NA      NA
#> --------------------------------------------
hessian(mlf)
#>           [,1] [,2]
#> [1,] -100.0018   NA
#> [2,]        NA   NA
# Note that now NA-s are in place of passive
# parameters.
# now invert only the free parameter part of the Hessian
sqrt(-solve(hessian(mlf)[activePar(mlf), activePar(mlf)]))
#>            [,1]
#> [1,] 0.09999911
# gives the standard deviation for the mean