// Laplace approximation from scratch demonstrated on 'spatial' example.
// The purpose is to show how the namespace 'autodiff' can be used to
// differentiate an iterative solver (the inner problem of the Laplace
// approximation). It is not meant to replace the approach used in the
// 'spatial' example (!)
/* Generic Laplace approximation (everywhere convex problems only). */
template<class Type, class Functor>
struct laplace_t {
Functor f; // User's implementation of joint nll
int niter; // Number of Newton iterations
f(f_), u(u_), niter(niter_) {}
Type operator()(){
// Solve inner problem - Newton iterations
for (int i=0; i<niter; i++){
}
// Laplace approximation
ans -= .5 * Type(u.size()) * log(2.0 * M_PI);
return ans;
}
};
template<class Type, class Functor>
laplace_t<Type, Functor> L(f, u, niter);
return L();
}
/* The following is (almost) copy-pasted from the 'spatial' example */
template<class Type>
struct joint_nll {
/* Data and parameter objects for spatial example: */
Type a;
Type log_sigma;
/* Constructor */
Type a_,
Type log_sigma_) :
y(y_), X(X_), dd(dd_), b(b_),
a(a_), log_sigma(log_sigma_) {}
/* Evaluate the negative joint log-likelihood as function of the
random effects */
template <typename T>
int n = u.size();
T res=0;
eta = eta + mu.template cast<T>();
for (int i=0; i<n; i++)
{
cov(i,i) = 1.0;
for (int j=0; j<i; j++)
{
// Exponentially decaying correlation
cov(i,j) = exp(-a * dd(i,j));
cov(j,i) = cov(i,j);
}
}
res += neg_log_density(u);
// logdpois = N log lam - lam
for(int i=0; i<y.size(); i++)
res -= T(y[i]) * eta[i] - exp(eta[i]);
return res;
}
};
template<class Type>
Type objective_function<Type>::operator() ()
{
int n = dd.rows();
// Construct joint negative log-likelihood
joint_nll<Type> jnll(y, X, dd, b, a, log_sigma);
// Random effect initial guess
u.setZero();
// Calculate Laplace approx (updates u)
Type res = laplace(jnll, u, niter);
return res;
}