multinomial_lpmf.hpp

Go to the documentation of this file.

#ifndef MULTINOMIAL_LPMF
#define MULTINOMIAL_LPMF
 
#include "density_components_base.hpp"
#include "../../common/fims_vector.hpp"
 
namespace fims_distributions {
template <typename Type>
struct MultinomialLPMF : public DensityComponentBase<Type> {
  fims::Vector<size_t> dims;

  MultinomialLPMF() : DensityComponentBase<Type>() {}

  virtual ~MultinomialLPMF() {}

  virtual const Type evaluate() {
    // set dims using data_observed_values if no user input
    if (dims.size() != 2) {
      dims.resize(2);
      dims[0] = this->data_observed_values->get_imax();
      dims[1] = this->data_observed_values->get_jmax();
    }
 
    // setup vector for recording the log probability density function values
    this->lpdf = static_cast<Type>(0.0); 
    this->lpdf_vec.resize(dims[0] * dims[1]);
    std::fill(this->lpdf_vec.begin(), this->lpdf_vec.end(), 0);
    size_t lpdf_vec_idx = 0; 
    // Dimension checks
    if (this->input_type == "data") {
      if (this->data_expected_values) {
        if (dims[0] * dims[1] != this->data_expected_values->size()) {
          throw std::invalid_argument(
              "MultinomialLPDF: Vector index out of bounds. The dimension of "
              "the "
              "number of rows times the number of columns is of size " +
              std::to_string(dims[0] * dims[1]) +
              " and the expected vector is of size " +
              std::to_string(this->data_expected_values->size()));
        }
      }
    } else {
      if (dims[0] * dims[1] != this->observed_values.size()) {
        throw std::invalid_argument(
            "MultinomialLPDF: Vector index out of bounds. The dimension of the "
            "number of  rows times the number of columns is of size " +
            std::to_string(dims[0] * dims[1]) +
            " and the observed vector is of size " +
            std::to_string(this->observed_values.size()));
      }
      if (this->observed_values.size() != this->expected_values.size()) {
        throw std::invalid_argument(
            "MultinomialLPDF: Vector index out of bounds. The dimension of the "
            "observed vector of size " +
            std::to_string(this->observed_values.size()) +
            " and the expected vector is of size " +
            std::to_string(this->expected_values.size()));
      }
    }
 
    for (size_t i = 0; i < dims[0]; i++) {
      // for each row, create new observed_values and prob vectors
      fims::Vector<Type> observed_values_vector;
      fims::Vector<Type> prob_vector;
      observed_values_vector.resize(dims[1]);
      prob_vector.resize(dims[1]);
 
      // Skips the entire row if any values are NA
      bool containsNA = false;
 
#ifdef TMB_MODEL
      for (size_t j = 0; j < dims[1]; j++) {
        if (this->input_type == "data") {
          // if data, check if there are any NA values and skip lpdf calculation
          // for entire row if there are
          if (this->get_observed(static_cast<size_t>(i),
                                 static_cast<size_t>(j)) ==
              this->data_observed_values->na_value) {
            containsNA = true;
            break;
          }
          if (!containsNA) {
            size_t idx = (i * dims[1]) + j;
            observed_values_vector[j] = this->get_observed(i, j);
            prob_vector[j] = this->get_expected(idx);
          }
        } else {
          // if not data (i.e. prior or process), use observed_values vector
          // instead of data_observed_values
          size_t idx = (i * dims[1]) + j;
          observed_values_vector[j] = this->get_observed(idx);
          prob_vector[j] = this->get_expected(idx);
        }
      }
 
      if (!containsNA) {
        Type row_lpmf = dmultinom(observed_values_vector.to_tmb(),
                                  prob_vector.to_tmb(), true);
        // The row log-probability mass has no unique split across bins because
        // the log(n!) term belongs to the whole row. Each bin reports an equal
        // share so that summing lpdf_vec over a composition counts each row
        // once.
        std::fill(this->lpdf_vec.begin() + lpdf_vec_idx,
                  this->lpdf_vec.begin() + lpdf_vec_idx + dims[1],
                  row_lpmf / static_cast<Type>(dims[1]));
 
        this->lpdf += row_lpmf;
      }
      lpdf_vec_idx += dims[1];
/*
if (this->simulate_flag)
{
    FIMS_SIMULATE_F(this->of)
    {
        fims::Vector<Type> sim_observed;
        sim_observed.resize(dims[1]);
        sim_observed = rmultinom(prob_vector);
        sim_observed.resize(this->observed_values.size());
        for (size_t j = 0; j < dims[1]; j++)
        {
            idx = (i * dims[1]) + j;
            this->observed_values[idx] = sim_observed[j];
        }
    }
}
*/
#endif
    }
 
#ifdef TMB_MODEL
#endif
    return (this->lpdf);
  }
};
}  // namespace fims_distributions
#endif