wishart.hpp

Go to the documentation of this file.

#ifndef __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__WISHART_HPP__
#define __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__WISHART_HPP__

#include <stan/prob/constants.hpp>
#include <stan/math/matrix_error_handling.hpp>
#include <stan/math/error_handling.hpp>
#include <stan/math/special_functions.hpp>
#include <stan/agrad/matrix.hpp>
#include <stan/prob/traits.hpp>
#include <boost/concept_check.hpp>

namespace stan {

  namespace prob {

    // Wishart(Sigma|n,Omega)  [Sigma, Omega symmetric, non-neg, definite; 
    //                          Sigma.dims() = Omega.dims();
    //                           n > Sigma.rows() - 1]
    template <bool propto,
              typename T_y, typename T_dof, typename T_scale, 
              class Policy>
    typename boost::math::tools::promote_args<T_y,T_dof,T_scale>::type
    wishart_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& W,
                const T_dof& nu,
                const Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic>& S,
                const Policy&) {
      static const char* function = "stan::prob::wishart_log(%1%)";

      using stan::math::check_greater_or_equal;
      using stan::math::check_size_match;
      using boost::math::tools::promote_args;

      typename Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic>::size_type k = W.rows();
      typename promote_args<T_y,T_dof,T_scale>::type lp(0.0);
      if (!check_greater_or_equal(function, nu, k - 1, 
                                  "Degrees of freedom parameter", &lp, Policy()))
        return lp;
      if (!check_size_match(function, 
                            W.rows(), "Rows of random variable",
                            W.cols(), "columns of random variable",
                            &lp, Policy()))
        return lp;
      if (!check_size_match(function, 
                            S.rows(), "Rows of scale parameter",
                            S.cols(), "columns of scale parameter",
                            &lp, Policy()))
        return lp;
      if (!check_size_match(function, 
                            W.rows(), "Rows of random variable",
                            S.rows(), "columns of scale parameter",
                            &lp, Policy()))
        return lp;
      // FIXME: domain checks

      Eigen::LLT< Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> > LLT_W = W.llt();
      if (LLT_W.info() != Eigen::Success) {
        lp = stan::math::policies::raise_domain_error<T_y>(function,
                                              "W is not positive definite (%1%)",
                                              0,Policy());
        return lp;
      }
      Eigen::LLT< Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic> > LLT_S = S.llt();
      if (LLT_S.info() != Eigen::Success) {
        lp = stan::math::policies::raise_domain_error<T_scale>(function,
                                              "S is not positive definite (%1%)",
                                              0,Policy());
        return lp;
      }

      Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> L_W = LLT_W.matrixL();
      Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic> L_S = LLT_S.matrixL();

      using stan::math::dot_product;
      using stan::math::mdivide_left_tri_low;
      using stan::math::crossprod;
      using stan::math::lmgamma;
      if (include_summand<propto,T_dof>::value)
        lp += nu * k * NEG_LOG_TWO_OVER_TWO;

      if (include_summand<propto,T_dof>::value)
        lp -= lmgamma(k, 0.5 * nu);

      if (include_summand<propto,T_dof,T_scale>::value)
        lp -= nu * L_S.diagonal().array().log().sum();

      if (include_summand<propto,T_scale,T_y>::value) {
        L_S = crossprod(mdivide_left_tri_low(L_S));
        Eigen::Matrix<T_scale,Eigen::Dynamic,1> S_inv_vec = Eigen::Map<
          const Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic> >(
                                &L_S(0), L_S.size(), 1);
        Eigen::Matrix<T_y,Eigen::Dynamic,1> W_vec = Eigen::Map<
          const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> >(
                                &W(0), W.size(), 1);
        lp -= 0.5 * dot_product(S_inv_vec, W_vec); // trace(S^-1 * W)
      }

      if (include_summand<propto,T_y,T_dof>::value && nu != (k + 1))
        lp += (nu - k - 1.0) * L_W.diagonal().array().log().sum();
      return lp;
    }

    template <bool propto,
              typename T_y, typename T_dof, typename T_scale>
    inline
    typename boost::math::tools::promote_args<T_y,T_dof,T_scale>::type
    wishart_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& W,
                const T_dof& nu,
                const Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic>& S) {
      return wishart_log<propto>(W,nu,S,stan::math::default_policy());
    }


    template <typename T_y, typename T_dof, typename T_scale, 
              class Policy>
    inline
    typename boost::math::tools::promote_args<T_y,T_dof,T_scale>::type
    wishart_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& W,
                const T_dof& nu,
                const Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic>& S,
                const Policy&) {
      return wishart_log<false>(W,nu,S,Policy());
    }


    template <typename T_y, typename T_dof, typename T_scale>
    inline
    typename boost::math::tools::promote_args<T_y,T_dof,T_scale>::type
    wishart_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& W,
                const T_dof& nu,
                const Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic>& S) {
      return wishart_log<false>(W,nu,S,stan::math::default_policy());
    }


  }
}
#endif