<div dir="ltr"><div><div><div>Hello all<br><br></div></div>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? <br>
<br></div><div>i've found functions on the following cite useful: <a href="http://qfc.fw.msu.edu/userfun.asp#r3">http://qfc.fw.msu.edu/userfun.asp#r3</a> . But are there other wasys of extracting the Hessian as well?<br>
</div><div><br></div>(Once one has the hessian in terms of the transformed parameters, it is eay to calculate the hessian in terms of the initial variables.)<br><div><br>regards<br><br><br><br><br><br>read.admbFit()<a name="r3"></a>
<pre style="font-family:Andale Mono,Lucida Console,Monaco,fixed,monospace;color:rgb(0,0,0);background-color:rgb(238,238,238);font-size:12px;border:1px dashed rgb(153,153,153);line-height:14px;padding:5px;overflow:auto;width:100%">
<code>
## 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)
}
</code></pre><br></div></div>