Maximum likelihood estimation

This is the main interface for the maxLik package, and the function that performs Maximum Likelihood estimation. It is a wrapper for different optimizers returning an object of class "maxLik". Corresponding methods handle the likelihood-specific properties of the estimates, including standard errors.

maxLik(logLik, grad = NULL, hess = NULL, start, method,
constraints=NULL, ...)

Arguments

logLik

log-likelihood function. Must have the parameter vector as the first argument. Must return either a single log-likelihood value, or a numeric vector where each component is log-likelihood of the corresponding individual observation.

grad

gradient of log-likelihood. Must have the parameter vector as the first argument. Must return either a single gradient vector with length equal to the number of parameters, or a matrix where each row is the gradient vector of the corresponding individual observation. If NULL, numeric gradient will be used.

hess

hessian of log-likelihood. Must have the parameter vector as the first argument. Must return a square matrix. If NULL, numeric Hessian will be used.

start

numeric vector, initial value of parameters. If it has names, these will also be used for naming the results.

method

maximisation method, currently either "NR" (for Newton-Raphson), "BFGS" (for Broyden-Fletcher-Goldfarb-Shanno), "BFGSR" (for the BFGS algorithm implemented in R), "BHHH" (for Berndt-Hall-Hall-Hausman), "SANN" (for Simulated ANNealing), "CG" (for Conjugate Gradients), or "NM" (for Nelder-Mead). Lower-case letters (such as "nr" for Newton-Raphson) are allowed. The default method is "NR" for unconstrained problems, and "NM" or "BFGS" for constrained problems, depending on if the grad argument was provided. "BHHH" is a good alternative given the likelihood is returned observation-wise (see maxBHHH).

Note that stochastic gradient ascent (SGA) is currently not supported as this method seems to be rarely used for maximum likelihood estimation.

constraints

either NULL for unconstrained maximization or a list, specifying the constraints. See maxBFGS.

...

further arguments, such as control, iterlim, or tol, are passed to the selected maximisation routine, i.e. maxNR, maxBFGS, maxBFGSR, maxBHHH, maxSANN, maxCG, or maxNM (depending on argument method). Arguments not used by the optimizers are forwarded to logLik, grad and hess.

Details

maxLik supports constrained optimization in the sense that constraints are passed further to the underlying optimization routines, and suitable default method is selected. However, no attempt is made to correct the resulting variance-covariance matrix. Hence the inference may be wrong. A corresponding warning is issued by the summary method.

Value

object of class 'maxLik' which inherits from class 'maxim'. Useful methods include

AIC: estimated parameter value
coef: estimated parameter value
logLik: log-likelihood value
nIter: number of iterations
stdEr: standard errors
summary: print summary table with estimates, standard errors, p, and z-values.
vcov: variance-covariance matrix

Warning

The constrained maximum likelihood estimation should be considered experimental. In particular, the variance-covariance matrix is not corrected for constrained parameter space.

Author

Ott Toomet, Arne Henningsen

Examples

## 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))
#> ----- Initial parameters: -----
#> fcn value: -45.47073 
#>      parameter initial gradient free
#> [1,]         1         54.52927    1
#> Condition number of the (active) hessian: 1 
#> -----Iteration 1 -----
#> -----Iteration 2 -----
#> -----Iteration 3 -----
#> -----Iteration 4 -----
#> -----Iteration 5 -----
#> --------------
#> gradient close to zero (gradtol) 
#> 5  iterations
#> estimate: 2.199217 
#> Function value: -21.18987 
summary( a )
#> --------------------------------------------
#> Maximum Likelihood estimation
#> Newton-Raphson maximisation, 5 iterations
#> Return code 1: gradient close to zero (gradtol)
#> Log-Likelihood: -21.18987 
#> 1  free parameters
#> Estimates:
#>      Estimate Std. error t value Pr(> t)    
#> [1,]    2.199      0.220   9.996  <2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> --------------------------------------------
##
## Estimate with analytic gradient and hessian.
## require much smaller tolerance
## setting 'tol=0' or negative essentially disables this stopping criterion
a <- maxLik(loglik, gradlik, hesslik, start=1,
            control=list(tol=-1, reltol=1e-12, gradtol=1e-12))
summary( a )
#> --------------------------------------------
#> Maximum Likelihood estimation
#> Newton-Raphson maximisation, 6 iterations
#> Return code 8: successive function values within relative tolerance limit (reltol)
#> Log-Likelihood: -21.18987 
#> 1  free parameters
#> Estimates:
#>      Estimate Std. error t value Pr(> t)    
#> [1,]   2.1992     0.2199      10  <2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> --------------------------------------------
##
## Next, we give an example with vector argument:
## fit normal distribution by estimating mean and standard deviation
## by maximum likelihood
##
loglik <- function(param) {
                           # param: vector of 2, c(mean, standard deviation)
   mu <- param[1]
   sigma <- param[2]
   ll <- -0.5*N*log(2*pi) - N*log(sigma) - sum(0.5*(x - mu)^2/sigma^2)
                           # can use dnorm(x, mu, sigma, log=TRUE) instead
   ll
}
x <- rnorm(100, 1, 2) # use mean=1, stdd=2
N <- length(x)
res <- maxLik(loglik, start=c(0,1)) # use 'wrong' start values
summary(res)
#> --------------------------------------------
#> Maximum Likelihood estimation
#> Newton-Raphson maximisation, 7 iterations
#> Return code 1: gradient close to zero (gradtol)
#> Log-Likelihood: -212.8427 
#> 2  free parameters
#> Estimates:
#>      Estimate Std. error t value  Pr(> t)    
#> [1,]   0.9775     0.2032    4.81 1.51e-06 ***
#> [2,]   2.0330     0.1438   14.14  < 2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> --------------------------------------------
##
## Same example, but now with named parameters and a fixed value
##
resFix <- maxLik(loglik, start=c(mu=0, sigma=1), fixed="sigma")
summary(resFix)  # 'sigma' is exactly 1.000 now.
#> --------------------------------------------
#> Maximum Likelihood estimation
#> Newton-Raphson maximisation, 2 iterations
#> Return code 8: successive function values within relative tolerance limit (reltol)
#> Log-Likelihood: -298.5382 
#> 1  free parameters
#> Estimates:
#>       Estimate Std. error t value Pr(> t)    
#> mu     0.97748    0.09998   9.777  <2e-16 ***
#> sigma  1.00000    0.00000      NA      NA    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> --------------------------------------------