diff --git a/TMB/inst/include/tmbutils/density.cpp b/TMB/inst/include/tmbutils/density.cpp index fa322f01d..d6596ddda 100644 --- a/TMB/inst/include/tmbutils/density.cpp +++ b/TMB/inst/include/tmbutils/density.cpp @@ -11,6 +11,7 @@ typedef matrix matrixtype; \ typedef array arraytype; #define VARIANCE_NOT_YET_IMPLEMENTED vectortype variance(){}; +#define JACOBIAN_NOT_YET_IMPLEMENTED arraytype jacobian(arraytype x){}; /** \brief Multivariate normal distribution with user supplied covariance matrix @@ -201,6 +202,59 @@ class N01{ VARIANCE_NOT_YET_IMPLEMENTED; }; + + +/** \brief Multivariate t distribution with user supplied scale matrix + + Class to evaluate the negative log density of a multivariate t distributed variable with general scale matrix Sigma and location vector 0 and df degrees of freedom. + This class should not be used as input distribution for SEPARABLE_t or PROJ_t. +*/ +template +class MVT_t: public MVNORM_t +{ + TYPEDEFS(scalartype_); + #include "lgamma.hpp" + scalartype df; + +public: + MVT_t() + : MVNORM_t() + {} + + MVT_t(scalartype df_) + : MVNORM_t() + { + setdf(df_); + } + MVT_t(matrixtype Sigma_, scalartype df_) + : MVNORM_t(Sigma_) + { + setdf(df_); + } + + void setdf(scalartype df_){ + df = df_; + } + + /** \brief Covariance extractor */ + matrixtype cov(){ + if(df > 2){ + return this->Sigma*df/(df-scalartype(2.0)); + } + } + + /** \brief Evaluate the negative log density */ + scalartype operator()(vectortype x){ + scalartype p = x.size(); + return -lgamma(scalartype(0.5)*(df+p))+lgamma(scalartype(0.5)*df)+p*scalartype(0.5)*log(df)+p*lgamma(scalartype(0.5))-scalartype(0.5)*this->logdetQ + scalartype(0.5)*(df+p)*log(scalartype(1.0)+this->Quadform(x)/df); + + } + JACOBIAN_NOT_YET_IMPLEMENTED; + VARIANCE_NOT_YET_IMPLEMENTED; +}; + + + /** \brief Stationary AR1 process Class to evaluate the negative log density of a (multivariate) AR1