# [ADMB Users] extract Hessian from within R

Ben Bolker bbolker at gmail.com
Mon Nov 4 06:27:01 PST 2013

```  Since this was cc'd to me:  I would recommend using read_pars() from
more) built in.  If you don't want to do that, you can go to

and find the code for reading the Hessian file:

{
nopar <- readBin(filen, what = "integer", n = 1)
hes <- readBin(filen, what = "double", n = nopar * nopar)
hes <- matrix(hes, byrow = TRUE, ncol = nopar)
colnames(hes) <- rownames(hes) <- sdparnames[seq(ncol(vcov))]
close(filen)
}

cheers
Ben Bolker

On 13-11-04 04:40 AM, Allan Clark wrote:
> 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
>
>
>
>
>
>
> |
> ## Function to read AD Model Builder fit from its par and cor files.
> ## Use by:
> ## Then the object 'simple.fit' is a list containing sub-objects
> # 'names', 'est', 'std', 'cor', and 'cov' for all model
> # parameters and sdreport quantities.
> #
>   ret<-list()
>   parfile<-as.numeric(scan(paste(file,'.par', sep=''),
>     what='', n=16, quiet=TRUE)[c(6,11,16)])
>   ret\$nopar<-as.integer(parfile)
>   ret\$nloglike<-parfile #objective function value
>
>   file<-paste(file,'.cor', sep='')
>   # total parameter including sdreport variables
>   ret\$totPar<-length(lin)-2
>   #log of the determinant of the hessian
>   ret\$logDetHess<-as.numeric(strsplit(lin, '=')[])
>   sublin<-lapply(strsplit(lin[1:ret\$totPar+2], ' '),function(x)x[x!=''])
>   ret\$names<-unlist(lapply(sublin,function(x)x))
>   ret\$est<-as.numeric(unlist(lapply(sublin,function(x)x)))
>   ret\$std<-as.numeric(unlist(lapply(sublin,function(x)x)))
>   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)
> }
> |
>
>

```