[ADMB Users] extract Hessian from within R

Cole Monnahan monnahc at uw.edu
Mon Nov 4 08:36:18 PST 2013


Allan,

There are two bounding functions controlled by the variable
"hybrid_bounded_flag", so if you want to know the transformation used for a
model you need to know that value. Fortunately it's included in the .hes
file. You can read more about this at:

http://www.admb-project.org/examples/admb-tricks/covariance-calculations

Hope that helps.

Cheers,
Cole



> Message: 2
> Date: Mon, 4 Nov 2013 11:40:47 +0200
> From: Allan Clark <allan.northpine at gmail.com>
> To: Admb users <Users at admb-project.org>, Ben Bolker
>         <bbolker at gmail.com>
> Subject: [ADMB Users]  extract Hessian from within R
> Message-ID:
>         <CAHpUJXDyNduANc0=5HTaywdaT03tNBUGaJ=
> kfzkD7iaFcnpW0g at mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
>
> 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-0001.html
> >
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.admb-project.org/pipermail/users/attachments/20131104/758a3845/attachment.html>


More information about the Users mailing list