search

kaskr / adcomp

AD computation with Template Model Builder (TMB)
Other
175 stars 80 forks source link

Delta-method for expectation of derived variables #26

Open James-Thorson opened 10 years ago

James-Thorson commented 10 years ago

Hi all,

It appears that sdreport() uses the delta-method to approximate the variance of derived values (specified as ADREPORT() outputs), but does not use the delta-method to approximate their expected values (and instead just reports f( E(theta) ) rather than E( f(theta) ), where theta are the estimated parameters, E(theta) are their MLE estimates, and f(theta) is the ADREPORT() computations).

Is it possible to add this computation? On first blush, it seems feasible given your AD calculation of the hessian -- but I don't fully understand the implication of the "generalized delta-method" computations used by sdreport().

cheers, Jim

skaug commented 10 years ago

Hi,

You are right that there is a conceptual issue when the model contains both (fixed effects) parameters theta and latent random variables "u". TMB calculates the MLE theta^ and the posterior mode estimates u^. For a function it returns f(theta^,u^) as the point estimate of f(theta,u). This is not a posterior mode estimate, but rather a hybrid approach.

It is worth thinking about if it is possible to return something that can be interpreted as a posterior mode.

Hans

James-Thorson commented 10 years ago

Hans,

I'm not sure I understand what "hybrid approach" means here, or how <theta^,mu^> differs from the posterior mode. I interpret that theta^ is the mode of the marginal likelihood, and mu^ is the mode of the conditional probability (i.e., uses Empirical Bayes), so please tell me if I'm wrong.

Lets say phi^ = <theta^,mu^>, and g^ = f(phi^), where g^ is a vector. if we can compute d2f/dphi2, it seems like we can calculate something similar to the delta-method approximation to E[g^]. I see how sdreport() uses:

obj2 <- MakeADFun(obj$env$data, obj$env$parameters, type="ADFun", ADreport=TRUE, DLL=obj$env$DLL) gr <- obj2$gr(par)

to calculate df/dphi. However, obj2$he seems to just return the hessian dNLL/dphi, where NLL is the objective function. Is there a way to configure TMB to calculate a hessian matrix for each ADREPORT element? I assume this would could be done individually for each element of g^ (i.e., calculating each E[g^] as being independent), and imagine this wouldn't require more memory than the original dNLL/dphi representation.

jim