TMB Documentation: multivariate_distributions.cpp

TMB Documentation v1.9.11

[フレーム]

multivariate_distributions.cpp

#include <TMB.hpp>

template<class Type>

Type objective_function<Type>::operator() ()

{

PARAMETER(dummy_par)

// Load namespace which contains the multivariate distributions

using namespace density;

Type res; // Dummy variable to which results will be asigned

// Example: MVNORM_t -----------------------------------------------------------

matrix<Type> 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> N_0_Sigma(Sigma); // N_0_Sigma is now a Distribution

res = N_0_Sigma(x0); // Evaluates (neg. log) density at x

res = MVNORM(Sigma)(x0); // Shorthand form of the above two lines

N_0_Sigma.cov(); // Returns covariance matrix (Sigma in this case)

REPORT(N_0_Sigma.cov()); // Report back to R

REPORT(MVNORM_t<Type>(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> nll(Lx);

vector<Type> x1(4);

x1.fill(2.0);

res = nll(x1 );

REPORT(nll.cov());

// 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> nllN01; // Distribution to be used at each time step

AR1_t<N01<Type> > nll2(phi,nllN01); // Create AR(1) process with N(0,1) distribution

REPORT(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(VECSCALE(nll2,sds2)(x2)); // Evaluate neg. log density

REPORT(VECSCALE(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> x(p,n); // Evaluation point

x.fill(1.0);

REPORT(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(AR1(phi,UNSTRUCTURED_CORR(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(AR1(phi)(x.transpose().col(0))); // First time series

REPORT(AR1(phi)(x.transpose().col(1))); // Second time series

return Type(0.0);

}

License: GPL v2