Open coraliewilliams opened 1 year ago
C++ code for multivariate distributions for different covariance structures:
#include <[TMB.hpp](http://kaskr.github.io/adcomp/TMB_8hpp.html)>
template<class Type>
Type objective_function<Type>::operator() ()
{
[PARAMETER](http://kaskr.github.io/adcomp/group__macros.html#ga70b2d449f7c53e87abee7d5fb71dc4c0)(dummy_par)
// Load namespace which contains the multivariate distributions
using namespace [density](http://kaskr.github.io/adcomp/namespacedensity.html);
Type res; // Dummy variable to which results will be asigned
// Example: MVNORM_t -----------------------------------------------------------
[matrix<Type>](http://kaskr.github.io/adcomp/structmatrix.html) Sigma(3,3);
Sigma.fill(0.1); // Fill the whole matrix
Sigma.diagonal() *= 10.0; // Multiply diagonal by 10 to positive definite Sigma
vector<Type> x0(3); // Point of evaluation
x0.fill(0.0); // Initialize x0 to be zero
[MVNORM_t<Type>](http://kaskr.github.io/adcomp/classdensity_1_1MVNORM__t.html) N_0_Sigma(Sigma); // N_0_Sigma is now a Distribution
res = N_0_Sigma(x0); // Evaluates (neg. log) density at x
res = [MVNORM](http://kaskr.github.io/adcomp/namespacedensity.html#a0721aacc998f4cbd0290626073c19763)(Sigma)(x0); // Shorthand form of the above two lines
N_0_Sigma.[cov](http://kaskr.github.io/adcomp/classdensity_1_1MVNORM__t.html#a1d284b3f3adac475399ae1e5428cb202)(); // Returns covariance matrix (Sigma in this case)
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)(N_0_Sigma.[cov](http://kaskr.github.io/adcomp/classdensity_1_1MVNORM__t.html#a1d284b3f3adac475399ae1e5428cb202)()); // Report back to R
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)([MVNORM_t<Type>](http://kaskr.github.io/adcomp/classdensity_1_1MVNORM__t.html)(Sigma).cov()); // Should return the same matrix as line above
// Example: UNSTRUCTURED_CORR_t -----------------------------------------------------
// First: read docs for UNSTRUCTURED_CORR_t
vector<Type> Lx(6); // Free parameters to be estimated
Lx.fill(2.5);
[UNSTRUCTURED_CORR_t<Type>](http://kaskr.github.io/adcomp/classdensity_1_1UNSTRUCTURED__CORR__t.html) nll(Lx);
vector<Type> x1(4);
x1.fill(2.0);
res = nll(x1 );
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)(nll.[cov](http://kaskr.github.io/adcomp/classdensity_1_1MVNORM__t.html#a1d284b3f3adac475399ae1e5428cb202)());
// Example: AR1_t -----------------------------------------------------
// Univariate case
int n = 10; // Number of time steps
vector<Type> x2(n); // Evaluation point, i.e. the time series
x2.fill(1.0);
Type phi = 0.5; // Autocorrelation
[N01<Type>](http://kaskr.github.io/adcomp/classdensity_1_1N01.html) nllN01; // Distribution to be used at each time step
[AR1_t<N01<Type>](http://kaskr.github.io/adcomp/classdensity_1_1AR1__t.html) > nll2(phi,nllN01); // Create AR(1) process with N(0,1) distribution
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)(nll2(x2)); // Evaluate neg. log density
AR1(phi)(x2); // Equivalent to line above
// Scale variances
vector<Type> sds2(n); // Vector of standard deviations (SD)
sds2.fill(2.0); // Set SD's
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)([VECSCALE](http://kaskr.github.io/adcomp/namespacedensity.html#ad4883b13b683e6b9e1e27846b74e4680)(nll2,sds2)(x2)); // Evaluate neg. log density
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)([VECSCALE](http://kaskr.github.io/adcomp/namespacedensity.html#ad4883b13b683e6b9e1e27846b74e4680)(AR1(phi),sds2)(x2)); // Equivalent to line above
// Multivariate AR(1) process where the innovation vector is correlated
int p=2; // dim(x)
[array<Type>](http://kaskr.github.io/adcomp/structarray.html) x(p,n); // Evaluation point
x.fill(1.0);
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)(x); // Reported back to R so that we can see layout of x
vector<Type> unconstrained_params(p*(p-1)/2);
unconstrained_params.fill(0.01); // Low correlaiton, but this is not directly the correlation
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)(AR1(phi,[UNSTRUCTURED_CORR](http://kaskr.github.io/adcomp/namespacedensity.html#aec3874ac89bc3a5d5d4d245db46b09e4)(unconstrained_params))(x)); // nll
// Find nll for univariate time series.
// Do not add to nll for system of time series due to intra series correlation.
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)(AR1(phi)(x.[transpose](http://kaskr.github.io/adcomp/structarray.html#a9f0df57cc5403f71194f2159ec80bd66)().col(0))); // First time series
[REPORT](http://kaskr.github.io/adcomp/group__macros.html#gac2d7888f4e0385930448f992de53079b)(AR1(phi)(x.[transpose](http://kaskr.github.io/adcomp/structarray.html#a9f0df57cc5403f71194f2159ec80bd66)().col(1))); // Second time series
return Type(0.0);
}
Reference to make matrix with Eigen in C++: http://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html
Links to documentation of covariance structures in glmmTMB: https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html http://kaskr.github.io/adcomp/classdensity_1_1UNSTRUCTURED__CORR__t.html