wishart_lpdf.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_PROB_WISHART_LPDF_HPP
#define STAN_MATH_PRIM_PROB_WISHART_LPDF_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/LDLT_factor.hpp>
#include <stan/math/prim/fun/lmgamma.hpp>
#include <stan/math/prim/fun/trace.hpp>
#include <stan/math/prim/fun/log_determinant_ldlt.hpp>
#include <stan/math/prim/fun/mdivide_left_ldlt.hpp>
#include <stan/math/prim/fun/constants.hpp>
 
namespace stan {
namespace math {
 
template <bool propto, typename T_y, typename T_dof, typename T_scale,
          require_stan_scalar_t<T_dof>* = nullptr,
          require_all_matrix_t<T_y, T_scale>* = nullptr>
return_type_t<T_y, T_dof, T_scale> wishart_lpdf(const T_y& W, const T_dof& nu,
                                                const T_scale& S) {
  using Eigen::Dynamic;
  using Eigen::Lower;
  using Eigen::Matrix;
  using T_W_ref = ref_type_t<T_y>;
  using T_nu_ref = ref_type_t<T_dof>;
  using T_S_ref = ref_type_t<T_scale>;
  static constexpr const char* function = "wishart_lpdf";
  Eigen::Index k = W.rows();
  check_size_match(function, "Rows of random variable", W.rows(),
                   "columns of scale parameter", S.rows());
 
  T_W_ref W_ref = W;
  T_nu_ref nu_ref = nu;
  T_S_ref S_ref = S;
 
  check_greater(function, "Degrees of freedom parameter", nu_ref, k - 1);
  check_square(function, "random variable", W_ref);
  check_square(function, "scale parameter", S_ref);
  check_symmetric(function, "random variable", W_ref);
  check_symmetric(function, "scale parameter", S_ref);
 
  auto ldlt_W = make_ldlt_factor(W_ref);
  check_ldlt_factor(function, "LDLT_Factor of random variable", ldlt_W);
  auto ldlt_S = make_ldlt_factor(S_ref);
  check_ldlt_factor(function, "LDLT_Factor of scale parameter", ldlt_S);
 
  return_type_t<T_y, T_dof, T_scale> lp(0.0);
 
  if (include_summand<propto, T_dof>::value) {
    lp -= nu_ref * k * HALF_LOG_TWO;
  }
 
  if (include_summand<propto, T_dof>::value) {
    lp -= lmgamma(k, 0.5 * nu_ref);
  }
 
  if (include_summand<propto, T_dof, T_scale>::value) {
    lp -= 0.5 * nu_ref * log_determinant_ldlt(ldlt_S);
  }
 
  if (include_summand<propto, T_scale, T_y>::value) {
    lp -= 0.5 * trace(mdivide_left_ldlt(ldlt_S, W_ref));
  }
 
  if (include_summand<propto, T_y, T_dof>::value && nu_ref != (k + 1)) {
    lp += 0.5 * (nu_ref - k - 1.0) * log_determinant_ldlt(ldlt_W);
  }
  return lp;
}
 
template <typename T_y, typename T_dof, typename T_scale>
inline return_type_t<T_y, T_dof, T_scale> wishart_lpdf(const T_y& W,
                                                       const T_dof& nu,
                                                       const T_scale& S) {
  return wishart_lpdf<false>(W, nu, S);
}
 
}  // namespace math
}  // namespace stan
#endif