[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