dirichlet.hpp

Go to the documentation of this file.

#ifndef __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__DIRICHLET_HPP__
#define __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__DIRICHLET_HPP__

#include <boost/math/special_functions/gamma.hpp>
#include <boost/random/gamma_distribution.hpp>
#include <boost/random/variate_generator.hpp>

#include <stan/prob/constants.hpp>
#include <stan/math/matrix_error_handling.hpp>
#include <stan/math/error_handling.hpp>
#include <stan/prob/traits.hpp>
#include <stan/math/functions/multiply_log.hpp>

namespace stan {

  namespace prob {

    template <bool propto,
              typename T_prob, typename T_prior_sample_size, 
              class Policy>
    typename boost::math::tools::promote_args<T_prob,T_prior_sample_size>::type
    dirichlet_log(const Eigen::Matrix<T_prob,Eigen::Dynamic,1>& theta,
              const Eigen::Matrix<T_prior_sample_size,Eigen::Dynamic,1>& alpha,
              const Policy&) {
      // FIXME: parameter check
      using boost::math::lgamma;
      using boost::math::tools::promote_args;
      typename promote_args<T_prob,T_prior_sample_size>::type lp(0.0);

      using stan::math::multiply_log;
      if (include_summand<propto,T_prior_sample_size>::value) {
        lp += lgamma(alpha.sum());
        for (int k = 0; k < alpha.rows(); ++k)
          lp -= lgamma(alpha[k]);
      }
      if (include_summand<propto,T_prob,T_prior_sample_size>::value)
        for (int k = 0; k < theta.rows(); ++k) 
          lp += multiply_log(alpha[k]-1, theta[k]);
      return lp;
    }

    template <bool propto,
              typename T_prob, typename T_prior_sample_size>
    inline
    typename boost::math::tools::promote_args<T_prob,T_prior_sample_size>::type
    dirichlet_log(const Eigen::Matrix<T_prob,Eigen::Dynamic,1>& theta,
          const Eigen::Matrix<T_prior_sample_size,Eigen::Dynamic,1>& alpha) {
      return dirichlet_log<propto>(theta,alpha,stan::math::default_policy());
    }


    template <typename T_prob, typename T_prior_sample_size, 
              class Policy>
    inline
    typename boost::math::tools::promote_args<T_prob,T_prior_sample_size>::type
    dirichlet_log(const Eigen::Matrix<T_prob,Eigen::Dynamic,1>& theta,
              const Eigen::Matrix<T_prior_sample_size,Eigen::Dynamic,1>& alpha,
              const Policy&) {
      return dirichlet_log<false>(theta,alpha,Policy());
    }

    template <typename T_prob, typename T_prior_sample_size>
    inline
    typename boost::math::tools::promote_args<T_prob,T_prior_sample_size>::type
    dirichlet_log(const Eigen::Matrix<T_prob,Eigen::Dynamic,1>& theta,
            const Eigen::Matrix<T_prior_sample_size,Eigen::Dynamic,1>& alpha) {
      return dirichlet_log<false>(theta,alpha,stan::math::default_policy());
    }

    template <class RNG>
    inline Eigen::VectorXd
    dirichlet_rng(const Eigen::Matrix<double,Eigen::Dynamic,1>& alpha,
                     RNG& rng) {
      using boost::variate_generator;
      using boost::gamma_distribution;

      double sum = 0;
      Eigen::VectorXd y(alpha.rows());
      for(int i = 0; i < alpha.rows(); i++) {
        variate_generator<RNG&, gamma_distribution<> >
          gamma_rng(rng, gamma_distribution<>(alpha(i,0),1));
        y(i) = gamma_rng();
        sum += y(i);
        }

      for(int i = 0; i < alpha.rows(); i++)
        y(i) /= sum;
      return y;
    }
  }
}
#endif