maxLik (version 0.6-0)

maxNR: Newton-Raphson Maximization

Description

Unconstrained maximization based on Newton-Raphson method.

Usage

maxNR(fn, grad = NULL, hess = NULL, start, print.level = 0,
      tol = 1e-08, reltol=sqrt(.Machine$double.eps), gradtol = 1e-06,
      steptol = 1e-10, lambdatol = 1e-06, qrtol = 1e-10,
      iterlim = 150,
      constraints=NULL,
      activePar=rep(TRUE, length(start)), ...)

Arguments

fn

function to be maximized. It must return either a single number or a numeric vector, which is summed (this is useful for BHHH method). If the parameters are out of range, fn should return NA. See details for consta

grad

gradient of the function. If NULL, numeric gradient is used. It must return a gradient vector, or matrix where columns correspond to individual parameters. The column sums are treated as gradient components (this is us

hess

Hessian matrix of the function. If missing, numeric Hessian, based on gradient, is used.

start

initial value for the parameter vector

print.level

this argument determines the level of printing which is done during the minimization process. The default value of 0 means that no printing occurs, a value of 1 means that initial and final details are printed and a value of 2 means that f

tol

stopping condition. Stop if the absolute difference between successive iterations less than tol, return code=2.

reltol

Relative convergence tolerance. The algorithm stops if it is unable to increase the value by a factor of 'reltol * (abs(val) + reltol)' at a step. Defaults to 'sqrt(.Machine$double.eps)', typically about '1e-8'.

gradtol

stopping condition. Stop if the norm of the gradient less than gradtol, return code=1.

steptol

stopping/error condition. If the quadratic approximation leads to lower function value instead of higher, or NA, the step length is halved and a new attempt is made. This procedure is repeated until step < steptol

lambdatol

control whether the Hessian is treated as negative definite. If the largest of the eigenvalues of the Hessian is larger than -lambdatol, a suitable diagonal matrix is subtracted from the Hessian (quadratic hill-climbing) in o

qrtol

QR-decomposition tolerance

iterlim

stopping condition. Stop if more than iterlim iterations, return code=4.

constraints

either NULL for unconstrained optimization or a list with two components eqA and eqB for equality-constrained optimization $A \theta + B = 0$. The constrained problem is forwarded to

activePar

index vector, which parameters are treated as free. Other parameters are treated as fixed constants.

...

further arguments to fn, grad and hess.

Value

• list of class "maxim" with following components:
• maximum
  fn value at maximum (the last calculated value if not converged).
• estimate
  estimated parameter value.
• gradient
  last gradient value which was calculated. Should be close to 0 if normal convergence.
• hessian
  Hessian at the maximum (the last calculated value if not converged).
• code
  return code:

  • 1

  {gradient close to zero (normal convergence).} 2{successive function values within tolerance limit (normal convergence).} 3{last step could not find higher value (probably not converged).} 4{iteration limit exceeded.} 5{Infinite value.} 6{Infinite gradient.} 7{Infinite Hessian.} 100{Initial value out of range.}

code

NULL

item

• message
• last.step
• f0
• climb
• activePar
• iterations
• type
• constraints
• outer.iterations
• barrier.value

itemize

• type

Details

The algorithm can work with constant parameters and related changes of parameter values. Constant parameters are useful if a parameter value is converging toward the boundary of support, or for testing.

One way is to put activePar to non-NULL, specifying which parameters we want to change. However, parameters can also be fixed in runtime by by signaling with fn. This may be useful if an estimation converges toward the edge of the parameter space possibly causing problems. The value of fn may have following attributes:

• constPar

{index vector. Which parameters are redefined to constant} newVal{a list with following components:

• index

{which parameters will have a new value} val{the new value of parameters} }

References

Greene, W., 2002, Advanced Econometrics Goldfeld, S. M. and Quandt, R. E., 1972, Nonlinear Methods in Econometrics. Amsterdam: North-Holland

See Also

nlm for Newton-Raphson optimization, optim for different gradient-based optimization methods.

Examples

## ML estimation of exponential duration model:
t <- rexp(100, 2)
loglik <- function(theta) sum(log(theta) - theta*t)
## Note the log-likelihood and gradient are summed over observations
gradlik <- function(theta) sum(1/theta - t)
hesslik <- function(theta) -100/theta^2
## Estimate with numeric gradient and Hessian
a <- maxNR(loglik, start=1, print.level=2)
summary(a)
## You would probably prefer 1/mean(t) instead ;-)
## Estimate with analytic gradient and Hessian
a <- maxNR(loglik, gradlik, hesslik, start=1)
summary(a)
##
## Next, we give an example with vector argument:  Estimate the mean and
## variance of a random normal sample by maximum likelihood
## Note: you might want to use maxLik instead
##
loglik <- function(param) {
  mu <- param[1]
  sigma <- param[2]
  ll <- -0.5*N*log(2*pi) - N*log(sigma) - sum(0.5*(x - mu)^2/sigma^2)
  ll
}
x <- rnorm(1000, 1, 2) # use mean=1, stdd=2
N <- length(x)
res <- maxNR(loglik, start=c(0,1)) # use 'wrong' start values
summary(res)
###
### Now an example of constrained optimization
###
## We maximize exp(-x^2 - y^2) where x+y = 1
f <- function(theta) {
  x <- theta[1]
  y <- theta[2]
  exp(-(x^2 + y^2))
  ## Note: you may want to use exp(- theta \%*\% theta) instead ;-)
}
## use constraints: x + y = 1
A <- matrix(c(1, 1), 1, 2)
B <- -1
res <- maxNR(f, start=c(0,0), constraints=list(eqA=A, eqB=B), print.level=1)
print(summary(res))