Nothing
##' Hessian (curvature matrix)
##'
##' Numerical evaluation of the Hessian of a real function f: \eqn{R^n }{R^n ->
##' R}\eqn{ \rightarrow R}{R^n -> R} on a generalized logit scale, i.e. using
##' transformed parameters according to x'=log((x-xl)/(xu-x))), with xl < x <
##' xu.
##'
##' This version uses a symmetric grid for the numerical evaluation computation
##' of first and second derivatives.
##'
##' @param x a list with components 'label' (of mode character), 'est' (the
##' parameter vector with the initial guess), 'low' (vector with lower bounds),
##' and 'upp' (vector with upper bounds)
##' @param f the function for which the Hessian is to be computed at point x
##' @param del step size on logit scale (numeric)
##' @param dapprox logical variable. If TRUE the off-diagonal elements are set
##' to zero. If FALSE (default) the full Hessian is computed
##' @param nfcn number of function calls
##' @return returns list with \item{df }{ first derivatives (logit scale) }
##' \item{ddf }{ Hessian (logit scale) } \item{nfcn }{ number of function calls
##' } \item{eigen }{ eigen values (logit scale) }
##' @note This function is part of the Bhat exploration tool
##' @author E. Georg Luebeck (FHCRC)
##' @seealso \code{\link{dfp}}, \code{\link{newton}}, \code{\link{ftrf}},
##' \code{\link{btrf}}, \code{\link{dqstep}}
##' @keywords array iteration methods optimize
##' @examples
##'
##' ## Rosenbrock Banana function
##' fr <- function(x) {
##' x1 <- x[1]
##' x2 <- x[2]
##' 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
##' }
##' ## define
##' x <- list(label=c("a","b"),est=c(1,1),low=c(-100,-100),upp=c(100,100))
##' logit.hessian(x,f=fr,del=dqstep(x,f=fr,sens=0.01))
##' ## shows the differences in curvature at the minimum of the Banana
##' ## function along principal axis (in a logit-transformed coordinate system)
##'
##' @export
##'
"logit.hessian" <-
function(x=x,f=f,del=rep(.002,length(x$est)),dapprox=FALSE,nfcn=0) {
# computes the hessian of f on logit scale
small <- 1.e-8
npar <- length(x$est)
nlm <- 2*npar*npar
np2 <- 2*npar
xt <- ftrf(x$est, x$low, x$upp)
f0 <- f(x$est); nfcn <- nfcn + 1
# cat('compute internal Hessian with function value at:',f0,'\n')
# *** INITIALIZE NEWTON BY CALCULATION OF A
# NLM POINT SIMPLEX GEOMETRY IN N DIM PARAMETER SPACE
xn <- matrix(rep(xt,nlm),ncol=nlm)
# *** ON AXES
# *** compute delta distance in each direction
for(i in 1:npar) {
j.even <- 2*i
j.odd <- j.even-1
xn[i,j.even] <- xn[i,j.even]-del[i]
xn[i,j.odd] <- xn[i,j.odd] +del[i]
}
if(dapprox==FALSE) {
# *** OFF AXIS
if(npar > 1) {
mc <- np2+1
for(i in 2:npar) {
for(j in 1:(i-1)) {
xn[i,mc] <- xn[i,mc]+del[i]
xn[j,mc] <- xn[j,mc]+del[j]
mc <- mc+1
xn[i,mc] <- xn[i,mc]-del[i]
xn[j,mc] <- xn[j,mc]+del[j]
mc <- mc+1
xn[i,mc] <- xn[i,mc]-del[i]
xn[j,mc] <- xn[j,mc]-del[j]
mc <- mc+1
xn[i,mc] <- xn[i,mc]+del[i]
xn[j,mc] <- xn[j,mc]-del[j]
mc <- mc+1
}}
}
f.vec <- numeric(nlm)
for(i in 1:nlm) {
f.vec[i] <- f(btrf(xn[,i], x$low, x$upp))
nfcn <- nfcn + 1
}
} else { # ignore off diagonal elements
f.vec <- numeric(np2)
for(i in 1:np2) {
f.vec[i] <- f(btrf(xn[,i], x$low, x$upp))
nfcn <- nfcn + 1
}
}
# *** FIRST AND DIAGONAL SECOND DERIVATIVES
i <- 1:npar; i.even <- 2*i; i.odd <- i.even-1
df <- (f.vec[i.even]-f.vec[i.odd])/2/del
if(npar > 1) {ddf <- diag((f.vec[i.even]+f.vec[i.odd]-2*f0)/(del**2))/2 }
else {ddf <- (f.vec[i.even]+f.vec[i.odd]-2*f0)/(del**2)/2 }
# print(format(ddf),quote=FALSE)
if(dapprox==FALSE) {
# *** SECOND DERIVATIVES
if(npar > 1) {
mc <- np2+1
for(i in 2:npar) {
for(j in 1:(i-1)) {
mc1 <- mc; mc2 <- mc1+1; mc3 <- mc2+1; mc4 <- mc3+1
ddf[i,j] <- (f.vec[mc1]+f.vec[mc3]-f.vec[mc2]-f.vec[mc4])/del[i]/del[j]/4
mc <- mc4+1
}}
ddf <- ddf+t(ddf)
}
} else {ddf <- 2 * ddf}
# *** EIGEN VALUES and MATRIX INVERSION
eig <- eigen(ddf)
# cat('eigen values: ',format(sort(eig$values)),'\n')
if(any(eig$values < 0)) {
warning('hessian not pos. definite')
}
if(any(abs(eig$values) < small)) {
warning('hessian may be singular')
}
# *** ADJUSTMENT OF del (feature not implemented)
# ddf.diag <- diag(ddf)
# del[ddf.diag > 0] <- .002/sqrt(ddf.diag[ddf.diag > 0])
return(list(df=df,ddf=ddf,nfcn=nfcn,eigen=eig$values))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.