Multivariate t distribution

TMB/inst/include/tmbutils/density.cpp

@@ -11,6 +11,7 @@ typedef matrix<scalartype> matrixtype;		\
 typedef array<scalartype> 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 scalartype_>
+class MVT_t: public MVNORM_t<scalartype_>
+{
+  TYPEDEFS(scalartype_);
+  #include "lgamma.hpp"
+  scalartype df;
+
+public:
+  MVT_t()
+    : MVNORM_t<scalartype>()
+  {}
+
+  MVT_t(scalartype df_)
+    : MVNORM_t<scalartype>()
+  {
+    setdf(df_);
+  }
+  MVT_t(matrixtype Sigma_, scalartype df_)
+    : MVNORM_t<scalartype>(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