Forum: help

RE: maxLik [ Reply ]
By: Arne Henningsen on 2016-03-22 21:53

1) maxLik() cannot be used "to find the hessian matrix"

2) In your example, the argument "x" of logLikHess(x) must be a vector of the (exactly two) parameters (not the x-values):
> logLikHess(c(1,2))
[,1] [,2]
[1,] -2.50 -5.2500
[2,] -5.25 -10.0625

3) A (less precise) numeric finite-difference calculation of the Hessian gives similar values:
> numericHessian( function(param) sum(dnorm(x,param[1],param[2],log=TRUE)), t0=c(1,2))
[,1] [,2]
[1,] -2.494005 -5.247358
[2,] -5.247358 -10.061285
which could indicate (although it does not proof) that the calculations in logLikHess() are correct.

RE: maxLik [ Reply ]
By: elebe nwezza on 2016-03-21 14:30

[forum:43068]

Please, am trying to understand how to use maxLik to find the hessian matrix of normal distribution. With the code below, is there need to call the maxLik function. If no, may you kindly help with example of how to use the maxLik

logLikHess <- function( param ) {
mu = param[ 1 ]
sigma = param[ 2 ]
N <- length( x )
logLikHessValues <- matrix( 0, nrow = 2, ncol = 2 )
logLikHessValues[ 1, 1 ] <- - N / sigma^2
logLikHessValues[ 1, 2 ] <- - 2 * sum( ( x - mu ) / sigma^3 )
logLikHessValues[ 2, 1 ] <- logLikHessValues[ 1, 2 ]
logLikHessValues[ 2, 2 ] <- N / sigma^2 - 3 * sum( ( x - mu )^2 / sigma^4 )
return(logLikHessValues)
}
x = c(1,2,4,2,3,1,5,3,5,5)
logLikHess(x)

Secondly, in the above code, is it correct to use x as the argument of logLikHess or a starting value such as c(2,1) as argument of logLikHess.

Thanks.

RE: maxLik [ Reply ] By: Arne Henningsen on 2016-03-19 10:09	[forum:43065]
As your model has two parameters (to be estimated), the function that returns the Hessian (in your case hesslik()) needs to return a symmetric 2x2 matrix. You can check this by calling this function, e.g. with the starting values: hesslik( c(3,1) ) If this does not return a symmetric 2x2 matrix, there is an error in hesslik().

RE: maxLik [ Reply ]
By: elebe nwezza on 2016-03-18 05:39

[forum:43062]

Thanks very much for your suggestion.
I tried following your suggestion with the code below but was flagged by an error message:

loglik = function(param) {
mu=param[1]
sigma=param[2]
x = c(3,2,4,1,5,4,6,8,7,9)
N = length(x)
ll = -0.5*N*log(2*pi) - N*log(sigma) - sum(0.5*(x - mu)^2/sigma^2)
return(ll)
}

res = maxLik2(loglik,gradlik, hesslik,start= c(3,1))

summary( res )

the error message is

Error in fn(start1, fixed = fixed, sumObs = TRUE, returnHessian = returnHessian, :
Wrong hessian dimension. Needed 2x2 but supplied 2x1

I tried the below

loglik = function(param) {
mu=param[1,1]
cv = param[1,2]
cva = cv
sigma=param[2,2]
x = c(3,2,4,1,5,4,6,8,7,9)
N = length(x)
ll = -0.5*N*log(2*pi) - N*log(sigma) - sum(0.5*(x - mu)^2/sigma^2)
return(ll)
}

res = maxLik2(loglik,gradlik, hesslik,start= c(3,1))

and the error message is

Error in param[1, 1] : incorrect number of dimensions

Pls, I need you help

RE: maxLik [ Reply ]
By: Arne Henningsen on 2016-03-17 12:13

[forum:43061]

Function maxLik() always requires argument "logLik," which must always be the log-likelihood function. In your command 'maxLik(hesslik, start=1)', argument "logLik" is specified as function hesslik(). However, function hesslik() (as defined in the code) does *not* return the log-likelihood value but the Hessian matrix. Therefore, your command 'maxLik(hesslik, start=1)' does not make sense. If you want to specify the Hessian matrix of the log-likelihood function, you can only provide it in *addition* to the log-likelihood function (and the gradient function):
a <- maxLik(loglik, gradlik, hesslik, start=1)

RE: maxLik [ Reply ]
By: elebe nwezza on 2016-03-17 09:10

[forum:43058]

Thanks for the example you gave me.

However, i observed that in computing the hessian which is suppose to be the inverse of variance of the estimator ( theta ) in the example you gave me, the result differ from the classical method unlike the mle.
though t is randomly generated

t <- rexp(100, 2)
loglik <- function(theta) log(theta) - theta*t
gradlik <- function(theta) 100/theta - t
hesslik <- function(theta) -100/theta^2

a <- maxLik2(loglik, start=1)
a
Maximum Likelihood estimation
Newton-Raphson maximisation, 5 iterations
Return code 1: gradient close to zero
Log-Likelihood: -26.27814 (1 free parameter(s))
Estimate(s): 2.090114

## the classical method
meant = solve(mean(t))
meant

[,1]
[1,] 2.090114

## for the hessian
t <- rexp(100, 2)
loglik <- function(theta) log(theta) - theta*t
gradlik <- function(theta) 100/theta - t
hesslik <- function(theta) -100/theta^2

a <- maxLik2(hesslik, start=1)

Maximum Likelihood estimation
Newton-Raphson maximisation, 150 iterations
Return code 4: Iteration limit exceeded.
Log-Likelihood: -0.0003644153 (1 free parameter(s))
Estimate(s): 523.8437

### the classical method
sigma = sd(t)^2
vart = length(t)/sigma
vart
[1] 414.8094

Pls, what may be the likely explanation for the difference in the hessian or my steps were wrong.

RE: maxLik [ Reply ]
By: Arne Henningsen on 2016-03-14 08:11

[forum:43044]

An example is given in the 'Example' section of maxLik's help-page:

## Estimate the parameter of exponential distribution
t <- rexp(100, 2)
loglik <- function(theta) log(theta) - theta*t
gradlik <- function(theta) 1/theta - t
hesslik <- function(theta) -100/theta^2

## Estimate with numeric gradient and hessian
a <- maxLik(loglik, start=1, control=list(printLevel=2))
summary( a )

## Estimate with analytic gradient and hessian
a <- maxLik(loglik, gradlik, hesslik, start=1)
summary( a )

maxLik [ Reply ] By: elebe nwezza on 2016-03-13 20:38	[forum:43043]
I want to provide the analytical derivative of the log-likelihood and then apply the maxLik function. Please, may so one help me out.