TMB Documentation  v1.9.1
TMBad/spde_epsilon.cpp
// Low-level demonstration of fast epsilon bias correction using 'sparse plus lowrank' hessian
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
using namespace R_inla;
using namespace density;
using namespace Eigen;
DATA_VECTOR(time);
DATA_IVECTOR(notcens);
DATA_IVECTOR(meshidxloc);
PARAMETER(log_tau);
PARAMETER(log_kappa);
PARAMETER(log_omega);
Type tau = exp(log_tau);
Type kappa = exp(log_kappa);
Type omega = exp(log_omega); // Parameter of Weibull distribution
Type nll = 0.0;
SparseMatrix<Type> Q = Q_spde(spde,kappa);
nll = GMRF(Q)(x); // Negative log likelihood
vector<Type> Xbeta = X*beta;
for(int i=0; i<time.size(); i++){
Type eta = Xbeta(i) + x(meshidxloc(i))/tau;
Type lambda = exp(eta);
Type t_omega = pow(time(i),omega);
Type S = exp(-lambda*t_omega); // Survival function
Type f = lambda*omega*t_omega/time(i)*S; // Weibull density
// Likelihood contribution depends on truncation status
if(notcens(i))
nll -= log(f);
else
nll -= log(S);
}
double nu = 1.0; // nu = alpha-d/2 = 2-1 by eqn (2) in Lindgren
Type rho = sqrt(8*nu)/kappa; // Distance at which correlation has dropped to 0.1 (p. 4 in Lindgren)
ADREPORT(rho);
// Optional: Epsilon contribution
PARAMETER_VECTOR(epsilon);
if (epsilon.size() > 0) {
Type S = x.exp().sum();
S = newton::Tag(S); // Set lowrank tag on S = sum(exp(x))
Type logS = log(S);
nll += epsilon[0] * logS;
}
return nll;
}
License: GPL v2