From 52aaec7a28f61d829f0142fea510e22cf07d5510 Mon Sep 17 00:00:00 2001 From: calbertsen Date: Mon, 8 Sep 2014 10:08:10 +0200 Subject: [PATCH 1/4] Added class to evaluate multivariate t distribution --- TMB/inst/include/tmbutils/density.cpp | 40 +++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/TMB/inst/include/tmbutils/density.cpp b/TMB/inst/include/tmbutils/density.cpp index fa322f01d..fd02701cd 100644 --- a/TMB/inst/include/tmbutils/density.cpp +++ b/TMB/inst/include/tmbutils/density.cpp @@ -201,6 +201,46 @@ 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. +*/ +template +class MVT_t: public MVNORM_t +{ + TYPEDEFS(scalartype_); + #include "lgamma.hpp" + scalartype df; + +public: + MVT_t(scalartype df_) + : MVNORM_t() + { + df = df_; + } + MVT_t(matrixtype Sigma_, scalartype df_) + : MVNORM_t(Sigma_) + { + df = df_; + } + + void setdf(scalartype df_){ + df = df_; + } + + /** \brief Evaluate the negative log density */ + scalartype operator()(vectortype x){ + scalartype p = x.size(); + //Lange et al. 1989 http://www.jstor.org/stable/2290063 + 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); + + } +}; + + + /** \brief Stationary AR1 process Class to evaluate the negative log density of a (multivariate) AR1 From b55b4b4741370fdb866951838952f585288115de Mon Sep 17 00:00:00 2001 From: calbertsen Date: Mon, 8 Sep 2014 10:19:02 +0200 Subject: [PATCH 2/4] Removed comment --- TMB/inst/include/tmbutils/density.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/TMB/inst/include/tmbutils/density.cpp b/TMB/inst/include/tmbutils/density.cpp index fd02701cd..57dd193e2 100644 --- a/TMB/inst/include/tmbutils/density.cpp +++ b/TMB/inst/include/tmbutils/density.cpp @@ -233,7 +233,6 @@ class MVT_t: public MVNORM_t /** \brief Evaluate the negative log density */ scalartype operator()(vectortype x){ scalartype p = x.size(); - //Lange et al. 1989 http://www.jstor.org/stable/2290063 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); } From f079a2fc47b009cdac19ae05e50142d32dd260d3 Mon Sep 17 00:00:00 2001 From: calbertsen Date: Mon, 8 Sep 2014 15:19:37 +0200 Subject: [PATCH 3/4] Corrected inherited functions --- TMB/inst/include/tmbutils/density.cpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/TMB/inst/include/tmbutils/density.cpp b/TMB/inst/include/tmbutils/density.cpp index 57dd193e2..4b5b25054 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(){}; /** \brief Multivariate normal distribution with user supplied covariance matrix @@ -205,7 +206,8 @@ class N01{ /** \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. + 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 other classes. */ template class MVT_t: public MVNORM_t @@ -218,24 +220,35 @@ class MVT_t: public MVNORM_t MVT_t(scalartype df_) : MVNORM_t() { - df = df_; + setdf(df_); } MVT_t(matrixtype Sigma_, scalartype df_) : MVNORM_t(Sigma_) { - df = df_; + setdf(df_); } void setdf(scalartype df_){ df = df_; } + /** \brief Covariance extractor */ + matrixtype cov(){ + if(df > 2){ + return Sigma*df/(df-scalartype(2.0)); + }else{ + return Sigma*scalartype(0.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; }; From a38a7f77c592d12891dae06dd35f7a5bb0cfc754 Mon Sep 17 00:00:00 2001 From: calbertsen Date: Tue, 9 Sep 2014 11:17:40 +0200 Subject: [PATCH 4/4] Added default constructor to MVT_t --- TMB/inst/include/tmbutils/density.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/TMB/inst/include/tmbutils/density.cpp b/TMB/inst/include/tmbutils/density.cpp index 4b5b25054..d6596ddda 100644 --- a/TMB/inst/include/tmbutils/density.cpp +++ b/TMB/inst/include/tmbutils/density.cpp @@ -11,7 +11,7 @@ typedef matrix matrixtype; \ typedef array arraytype; #define VARIANCE_NOT_YET_IMPLEMENTED vectortype variance(){}; -#define JACOBIAN_NOT_YET_IMPLEMENTED arraytype jacobian(){}; +#define JACOBIAN_NOT_YET_IMPLEMENTED arraytype jacobian(arraytype x){}; /** \brief Multivariate normal distribution with user supplied covariance matrix @@ -207,7 +207,7 @@ class N01{ /** \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 other classes. + This class should not be used as input distribution for SEPARABLE_t or PROJ_t. */ template class MVT_t: public MVNORM_t @@ -217,6 +217,10 @@ class MVT_t: public MVNORM_t scalartype df; public: + MVT_t() + : MVNORM_t() + {} + MVT_t(scalartype df_) : MVNORM_t() { @@ -235,10 +239,8 @@ class MVT_t: public MVNORM_t /** \brief Covariance extractor */ matrixtype cov(){ if(df > 2){ - return Sigma*df/(df-scalartype(2.0)); - }else{ - return Sigma*scalartype(0.0); - } + return this->Sigma*df/(df-scalartype(2.0)); + } } /** \brief Evaluate the negative log density */