[ADMB Users] extract Hessian from within R

Allan Clark allan.northpine at gmail.com
Mon Nov 4 01:40:47 PST 2013


Hello all

I've coded up a simple model where I estimate various probabilities. These
are restricted to (0,1) and decided to rather transform these variables so
that they are not bounded anymore. How does one extract the hessian of the
parameters from within R?

i've found functions on the following cite useful:
http://qfc.fw.msu.edu/userfun.asp#r3 . But are there other wasys of
extracting the Hessian as well?

(Once one has the hessian in terms of the transformed parameters, it is eay
to calculate the hessian in terms of the initial variables.)

regards





read.admbFit()


## Function to read AD Model Builder fit from its par and cor files.
## Use by:
## simple.fit <- read.admbFit('c:/admb/examples/simple')
## Then the object 'simple.fit' is a list containing sub-objects
# 'names', 'est', 'std', 'cor', and 'cov' for all model
# parameters and sdreport quantities.
#
read.admbFit<-function(file){
  ret<-list()
  # read par file
  parfile<-as.numeric(scan(paste(file,'.par', sep=''),
    what='', n=16, quiet=TRUE)[c(6,11,16)])
  ret$nopar<-as.integer(parfile[1])
  ret$nloglike<-parfile[2] #objective function value
  ret$maxgrad<-parfile[3]

  # read cor file
  file<-paste(file,'.cor', sep='')
  lin<-readLines(file)
  # total parameter including sdreport variables
  ret$totPar<-length(lin)-2
  #log of the determinant of the hessian
  ret$logDetHess<-as.numeric(strsplit(lin[1], '=')[[1]][2])
  sublin<-lapply(strsplit(lin[1:ret$totPar+2], ' '),function(x)x[x!=''])
  ret$names<-unlist(lapply(sublin,function(x)x[2]))
  ret$est<-as.numeric(unlist(lapply(sublin,function(x)x[3])))
  ret$std<-as.numeric(unlist(lapply(sublin,function(x)x[4])))
  ret$cor<-matrix(NA, ret$totPar, ret$totPar)
  corvec<-unlist(sapply(1:length(sublin), function(i)sublin[[i]][5:(4+i)]))
  ret$cor[upper.tri(ret$cor, diag=TRUE)]<-as.numeric(corvec)
  ret$cor[lower.tri(ret$cor)] <- t(ret$cor)[lower.tri(ret$cor)]
  # covariance matrix
  ret$cov<-ret$cor*(ret$std %o% ret$std)
  return(ret)
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.admb-project.org/pipermail/users/attachments/20131104/9cb46bc9/attachment.html>


More information about the Users mailing list