ordered_probit_lpmf.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_PROB_ORDERED_PROBIT_LPMF_HPP
#define STAN_MATH_PRIM_PROB_ORDERED_PROBIT_LPMF_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/log_diff_exp.hpp>
#include <stan/math/prim/fun/log1m_exp.hpp>
#include <stan/math/prim/fun/scalar_seq_view.hpp>
#include <stan/math/prim/fun/size_mvt.hpp>
#include <stan/math/prim/fun/to_ref.hpp>
#include <stan/math/prim/prob/std_normal_lcdf.hpp>
#include <stan/math/prim/fun/vector_seq_view.hpp>
#include <vector>
#include <cmath>
 
namespace stan {
namespace math {
 
template <bool propto, typename T_y, typename T_loc, typename T_cut>
return_type_t<T_loc, T_cut> ordered_probit_lpmf(const T_y& y,
                                                const T_loc& lambda,
                                                const T_cut& c) {
  using std::exp;
  using std::log;
  using T_lambda_ref = ref_type_t<T_loc>;
  static constexpr const char* function = "ordered_probit";
 
  check_nonzero_size(function, "Cut-points", c);
  int N = math::size(lambda);
  int C_l = size_mvt(c);
  vector_seq_view<T_cut> c_vec(c);
  int K = c_vec[0].size() + 1;
 
  check_consistent_sizes(function, "Integers", y, "Locations", lambda);
  if (C_l > 1) {
    check_size_match(function, "Length of location variables ", N,
                     "Number of cutpoint vectors ", C_l);
  }
 
  check_bounded(function, "Random variable", y, 1, K);
  scalar_seq_view<T_y> y_vec(y);
  check_nonzero_size(function, "First cutpoint set", c_vec[0]);
  for (int i = 0; i < size_mvt(c); ++i) {
    check_size_match(function, "One cutpoint set", K - 1, "First cutpoint set",
                     c_vec[i].size());
    check_ordered(function, "Cut-points", c_vec[i]);
    if (K > 2) {
      check_finite(function, "Final cut point", c_vec[i].coeff(K - 2));
    }
    check_finite(function, "First cut point", c_vec[i].coeff(0));
  }
 
  T_lambda_ref lambda_ref = lambda;
  check_finite(function, "Location parameter", lambda_ref);
  scalar_seq_view<T_lambda_ref> lambda_vec(lambda_ref);
 
  return_type_t<T_loc, T_cut> logp_n(0.0);
 
  for (int i = 0; i < N; ++i) {
    int K = c_vec[i].size() + 1;
 
    if (y_vec[i] == 1) {
      logp_n += std_normal_lcdf(c_vec[i].coeff(0) - lambda_vec[i]);
    } else if (y_vec[i] == K) {
      logp_n += std_normal_lcdf(lambda_vec[i] - c_vec[i].coeff(K - 2));
    } else {
      logp_n += log_diff_exp(
          std_normal_lcdf(c_vec[i].coeff(y_vec[i] - 1) - lambda_vec[i]),
          std_normal_lcdf(c_vec[i].coeff(y_vec[i] - 2) - lambda_vec[i]));
    }
  }
  return logp_n;
}
 
template <typename T_y, typename T_loc, typename T_cut>
return_type_t<T_loc, T_cut> ordered_probit_lpmf(const T_y& y,
                                                const T_loc& lambda,
                                                const T_cut& c) {
  return ordered_probit_lpmf<false>(y, lambda, c);
}
}  // namespace math
}  // namespace stan
#endif