score function

[Hervé DAKPO]

Dear all, I am doing maximum likelihood estimation over a sample of firms (using skew normal distribution). For my paper the Reviewers ask to compute robust standard error à la sandwich. This requires to obtain the gradient or score function for each observation (this is also necessary if one wishes to obtain BHHH hessian approximation).

The thing is the log-likelihood is the sum over observations. Does anybody know a way to obtain the score function for each observation?

Thanks

[Kasper Kristensen]

Based on your other post I assume this a RTMB question (not TMB).

An easy way to get gradients of each term is to ADREPORT all terms. Here, using a simple normal model:

## Example model and likelihood
library(RTMB)
nobs <- 10
parms <- list(mu=0, sd=1)
xobs <- rnorm(nobs, mean=2, sd=3)
f <- function(parms) {
    getAll(parms)
    nll <- -dnorm(xobs, mu, sd, log=TRUE)
    ADREPORT(nll)
    sum(nll)
}
## Get MLE
obj <- MakeADFun(f, parms)
opt <- nlminb(obj$par, obj$fn, obj$gr)

You can now get the gradient of each term (row) wrt to each parameter (column)like this:

obj <- MakeADFun(f, parms, ADreport=TRUE)
obj$gr(opt$par)

However, note that this approach is very slow when the number of terms (nobs) is much higher than the number of parameters (2) because TMB uses reverse mode AD.

If this is an issue, you can get the same gradients another way (using an 'adjoint trick'):

## From T: R^2 -> R^nobs

## Generate T2: R^nobs -> R^2

T <- GetTape(obj)
T2 <- MakeTape(function(weight) {
    WT <- MakeTape(function(x) sum(T(x) * weight), opt$par)
    WT$jacfun()(advector(opt$par))
}, rep(1,nobs))
t(T2$jacobian(rep(1,nobs)))

[James Thorson]

Kasper,

Cool trick, and thanks for sharing that!  Is there any way to do the same thing in TMB (rather than RTMB)?  

Obviously the same questions arise for both implementations, and it seems like useful machinery to repurpose for other uses too.

Jim

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/fcc573c3-5881-4ecd-8a07-69845a1d8426n%40googlegroups.com.

[Kasper Kristensen]

Jim, you can apply 'RTMB::GetTape(obj)' for a plain TMB model as well, so the code example should be generally applicable (using RTMB::MakeTape, RTMB::advector etc). Only restriction is that the TMB model must be binary compatible, i.e compiled with identical config variables/flags (e.g. AD framework, OpenMP, etc) as RTMB.

[James Thorson]

Wow, mind blown that we could use RTMB utilities for a TMB model :0

I ran `RTMB::GetTape(obj)` on a phylosem object (which already uses framework="TMBad"), and got a helpful message about updating TMB and RTMB.  I then installed the github TMB, the CRAN RTMB version, and reinstalled phylosem.  Trying `RTMB::GetTape(obj)` again, I get the errors below:

I'm of course happy to share a simple demo script.  

Most immediately, I'm excited to see if I can use the igrah plotting for the tape :0  But it seems like there's a ton of interesting uses to explore with this type of functionality!

To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/138d5eaa-ccac-4975-afff-b790dbb6da4cn%40googlegroups.com.

[Hervé DAKPO]

Thanks a lot the trick works perfectly.

Merci