TMB Documentation  v1.9.0
spde_aniso_speedup.cpp
// Speedup "spde_aniso.cpp" by moving normalization out of the template.
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
using namespace R_inla;
using namespace density;
using namespace Eigen;
DATA_INTEGER(flag); // flag=0 => only prior
DATA_VECTOR(time);
DATA_IVECTOR(notcens);
DATA_IVECTOR(meshidxloc);
PARAMETER(log_tau);
PARAMETER(log_kappa);
PARAMETER_VECTOR(ln_H_input);
PARAMETER(log_alpha);
Type tau = exp(log_tau);
Type kappa = exp(log_kappa);
Type alpha = exp(log_alpha);
Type nll = 0.0;
// Need to parameterize H matrix such that det(H)=1 (preserving volume)
// Note that H appears in (20) in Lindgren et al 2011
matrix<Type> H(2,2);
H(0,0) = exp(ln_H_input(0));
H(1,0) = ln_H_input(1);
H(0,1) = ln_H_input(1);
H(1,1) = (1+ln_H_input(1)*ln_H_input(1)) / exp(ln_H_input(0));
SparseMatrix<Type> Q = Q_spde(spde,kappa,H);
REPORT(H)
// Drop normalizing constant for random effects
nll = GMRF(Q, false)(x); // Negative log likelihood
// Return un-normalized density on request
if (flag == 0) return nll;
// Weibull likelihood with cencoring
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_alpha = pow(time(i),alpha);
Type S = exp(-lambda*t_alpha); // Survival function
Type f = lambda*alpha*t_alpha/time(i)*S; // Densities
// Likelihood contribution depends on truncation status
if(notcens(i))
nll -= log(f);
else
nll -= log(S);
}
return nll;
}
License: GPL v2