Beta Hierarchy

docs/protos.html

@@ -457,10 +457,6 @@ <h2>Table of Contents</h2>
                   <a href="#bayesmix.TruncSBPrior.DPPrior"><span class="badge">M</span>TruncSBPrior.DPPrior</a>
                 </li>
 
-                <li>
-                  <a href="#bayesmix.TruncSBPrior.MFMPrior"><span class="badge">M</span>TruncSBPrior.MFMPrior</a>
-                </li>
-
                 <li>
                   <a href="#bayesmix.TruncSBPrior.PYPrior"><span class="badge">M</span>TruncSBPrior.PYPrior</a>
                 </li>
@@ -2369,24 +2365,18 @@ <h3 id="bayesmix.TruncSBPrior">TruncSBPrior</h3>
                 </tr>
 
                 <tr>
-                  <td>mfm_prior</td>
-                  <td><a href="#bayesmix.TruncSBPrior.MFMPrior">TruncSBPrior.MFMPrior</a></td>
-                  <td></td>
-                  <td><p> </p></td>
-                </tr>
-
-                <tr>
-                  <td>pm_prior</td>
-                  <td><a href="#bayesmix.PythonMixPrior">PythonMixPrior</a></td>
+                  <td>num_components</td>
+                  <td><a href="#uint32">uint32</a></td>
                   <td></td>
-                  <td><p> </p></td>
+                  <td><p>Number of components in the process </p></td>
                 </tr>
 
                 <tr>
-                  <td>num_components</td>
-                  <td><a href="#uint32">uint32</a></td>
+                  <td>infinite_mixture</td>
+                  <td><a href="#bool">bool</a></td>
                   <td></td>
-                  <td><p>Number of components in the process </p></td>
+                  <td><p>If true we must use the Slice Sampler, and num_components is used only for
+the initialization </p></td>
                 </tr>
 
             </tbody>
@@ -2444,30 +2434,6 @@ <h3 id="bayesmix.TruncSBPrior.DPPrior">TruncSBPrior.DPPrior</h3>
 
 
 
-        <h3 id="bayesmix.TruncSBPrior.MFMPrior">TruncSBPrior.MFMPrior</h3>
-        <p></p>
-
-
-          <table class="field-table">
-            <thead>
-              <tr><td>Field</td><td>Type</td><td>Label</td><td>Description</td></tr>
-            </thead>
-            <tbody>
-
-                <tr>
-                  <td>totalmass</td>
-                  <td><a href="#double">double</a></td>
-                  <td></td>
-                  <td><p>Truncated Dirichlet process </p></td>
-                </tr>
-
-            </tbody>
-          </table>
-
-
-
-
-
         <h3 id="bayesmix.TruncSBPrior.PYPrior">TruncSBPrior.PYPrior</h3>
         <p></p>
 

src/hierarchies/CMakeLists.txt

@@ -8,6 +8,7 @@ target_sources(bayesmix
     lin_reg_uni_hierarchy.h
     fa_hierarchy.h
     lapnig_hierarchy.h
+    betagg_hierarchy.h
 )
 
 add_subdirectory(likelihoods)

src/hierarchies/betagg_hierarchy.h

@@ -0,0 +1,57 @@
+#ifndef BAYESMIX_HIERARCHIES_BETA_GG_HIERARCHY_H_
+#define BAYESMIX_HIERARCHIES_BETA_GG_HIERARCHY_H_
+
+#include "base_hierarchy.h"
+#include "hierarchy_id.pb.h"
+#include "likelihoods/beta_likelihood.h"
+#include "priors/gamma_gamma_prior.h"
+#include "updaters/random_walk_updater.h"
+
+/**
+ * Beta Gamma-Gamma hierarchy for univaraite data in [0, 1]
+ *
+ * This class represents a hierarchical model where data are distributed
+ * according to a Beta likelihood (see the `BetaLikelihood` class for
+ * details). The shape and rate parameters of the likelihood have
+ * independent gamma priors. That is
+ *
+ * \f[
+ *    f(x_i \mid \alpha, \beta) &= Beta(\alpha, \beta) \\
+ *    \alpha &\sim Gamma(\alpha_a, \alpha_b) \\
+ *    \beta &\sim Gamma(\beta_a, \beta_b)
+ * \f]
+ *
+ * The state is composed of shape and rate. Note that this hierarchy
+ * is NOT conjugate, meaning that the marginal distribution is not available
+ * in closed form.
+ */
+
+class BetaGGHierarchy
+    : public BaseHierarchy<BetaGGHierarchy, BetaLikelihood, GGPriorModel> {
+ public:
+  BetaGGHierarchy() = default;
+  ~BetaGGHierarchy() = default;
+
+  //! Returns the Protobuf ID associated to this class
+  bayesmix::HierarchyId get_id() const override {
+    return bayesmix::HierarchyId::BetaGG;
+  }
+
+  //! Sets the default updater algorithm for this hierarchy
+  void set_default_updater() {
+    updater = std::make_shared<RandomWalkUpdater>(0.1);
+  }
+
+  //! Initializes state parameters to appropriate values
+  void initialize_state() override {
+    // Get hypers
+    auto hypers = prior->get_hypers();
+    // Initialize likelihood state
+    State::ShapeRate state;
+    state.shape = hypers.a_shape / hypers.a_rate;
+    state.rate = hypers.b_shape / hypers.b_rate;
+    like->set_state(state);
+  };
+};
+
+#endif  // BAYESMIX_HIERARCHIES_BETA_GG_HIERARCHY_H_

src/hierarchies/likelihoods/CMakeLists.txt

@@ -12,6 +12,8 @@ target_sources(bayesmix PUBLIC
     laplace_likelihood.cc
     fa_likelihood.h
     fa_likelihood.cc
+    beta_likelihood.h
+    beta_likelihood.cc
 )
 
 add_subdirectory(states)

src/hierarchies/likelihoods/beta_likelihood.cc

@@ -0,0 +1,22 @@
+#include "beta_likelihood.h"
+
+double BetaLikelihood::compute_lpdf(const Eigen::RowVectorXd &datum) const {
+  return stan::math::beta_lpdf(datum(0), state.shape, state.rate);
+}
+
+void BetaLikelihood::update_sum_stats(const Eigen::RowVectorXd &datum,
+                                      bool add) {
+  double x = datum(0);
+  if (add) {
+    sum_logs += std::log(x);
+    sum_logs1m += std::log(1. - x);
+  } else {
+    sum_logs -= std::log(x);
+    sum_logs1m -= std::log(1. - x);
+  }
+}
+
+void BetaLikelihood::clear_summary_statistics() {
+  sum_logs = 0;
+  sum_logs1m = 0;
+}

src/hierarchies/likelihoods/beta_likelihood.h

@@ -0,0 +1,56 @@
+#ifndef BAYESMIX_HIERARCHIES_LIKELIHOODS_BETA_LIKELIHOOD_H_
+#define BAYESMIX_HIERARCHIES_LIKELIHOODS_BETA_LIKELIHOOD_H_
+
+#include <google/protobuf/stubs/casts.h>
+
+#include <memory>
+#include <stan/math/rev.hpp>
+#include <vector>
+
+#include "algorithm_state.pb.h"
+#include "base_likelihood.h"
+#include "states/includes.h"
+
+/**
+ * A univariate Beta likelihood, using the `State::ShapeRate` state. Represents
+ * the model:
+ *
+ * \f[
+ *    y_1,\dots,y_k \mid \mu, \sigma^2 \stackrel{\small\mathrm{iid}}{\sim}
+ * Beta(a, b),
+ * \f]
+ */
+
+class BetaLikelihood
+    : public BaseLikelihood<BetaLikelihood, State::ShapeRate> {
+ public:
+  BetaLikelihood() = default;
+  ~BetaLikelihood() = default;
+  bool is_multivariate() const override { return false; };
+  bool is_dependent() const override { return false; };
+  void clear_summary_statistics() override;
+
+  template <typename T>
+  T cluster_lpdf_from_unconstrained(
+      const Eigen::Matrix<T, Eigen::Dynamic, 1> &unconstrained_params) const {
+    assert(unconstrained_params.size() == 2);
+
+    T a = stan::math::positive_constrain(unconstrained_params(0));
+    T b = stan::math::positive_constrain(unconstrained_params(1));
+
+    T out = 0.;
+
+    return (a - 1.) * sum_logs + (b - 1.) * sum_logs1m -
+           card * (stan::math::lgamma(a) + stan::math::lgamma(b) -
+                   stan::math::lgamma(a + b));
+  }
+
+ protected:
+  double compute_lpdf(const Eigen::RowVectorXd &datum) const override;
+  void update_sum_stats(const Eigen::RowVectorXd &datum, bool add) override;
+
+  double sum_logs = 0;
+  double sum_logs1m = 0;
+};
+
+#endif  // BAYESMIX_HIERARCHIES_LIKELIHOODS_BETA_LIKELIHOOD_H_

src/hierarchies/likelihoods/states/includes.h

@@ -3,6 +3,7 @@
 
 #include "fa_state.h"
 #include "multi_ls_state.h"
+#include "shape_rate_state.h"
 #include "uni_lin_reg_ls_state.h"
 #include "uni_ls_state.h"
 

src/hierarchies/likelihoods/states/shape_rate_state.h

@@ -0,0 +1,88 @@
+#ifndef BAYESMIX_HIERARCHIES_LIKELIHOODS_STATES_SHAPE_RATE_STATE_H_
+#define BAYESMIX_HIERARCHIES_LIKELIHOODS_STATES_SHAPE_RATE_STATE_H_
+
+#include <memory>
+#include <stan/math/rev.hpp>
+
+#include "algorithm_state.pb.h"
+#include "base_state.h"
+#include "src/utils/proto_utils.h"
+
+namespace State {
+
+//! Returns the constrained parametrization from the
+//! unconstrained one, i.e. [in[0], exp(in[1])]
+template <typename T>
+Eigen::Matrix<T, Eigen::Dynamic, 1> shape_rate_to_constrained(
+    Eigen::Matrix<T, Eigen::Dynamic, 1> in) {
+  Eigen::Matrix<T, Eigen::Dynamic, 1> out(2);
+  out << stan::math::exp(in(0)), stan::math::exp(in(1));
+  return out;
+}
+
+//! Returns the unconstrained parametrization from the
+//! constrained one, i.e. [log(in[0]), log(in[1])]
+template <typename T>
+Eigen::Matrix<T, Eigen::Dynamic, 1> shape_rate_to_unconstrained(
+    Eigen::Matrix<T, Eigen::Dynamic, 1> in) {
+  Eigen::Matrix<T, Eigen::Dynamic, 1> out(2);
+  out << stan::math::log(in(0)), stan::math::log(in(1));
+  return out;
+}
+
+//! Returns the log determinant of the jacobian of the map
+//! (x, y) -> (log(x), log(y)), that is the inverse map of the
+//! constrained -> unconstrained representation.
+template <typename T>
+T shape_rate_log_det_jac(Eigen::Matrix<T, Eigen::Dynamic, 1> constrained) {
+  T out = 0;
+  stan::math::positive_constrain(stan::math::log(constrained(0)), out);
+  stan::math::positive_constrain(stan::math::log(constrained(1)), out);
+  return out;
+}
+
+//! A univariate shape-rate state
+//! The unconstrained representation corresponds to (log(shape), log(rate))
+class ShapeRate : public BaseState {
+ public:
+  double shape, rate;
+
+  using ProtoState = bayesmix::AlgorithmState::ClusterState;
+
+  Eigen::VectorXd get_unconstrained() const override {
+    Eigen::VectorXd temp(2);
+    temp << shape, rate;
+    return shape_rate_to_unconstrained(temp);
+  }
+
+  void set_from_unconstrained(const Eigen::VectorXd &in) override {
+    Eigen::VectorXd temp = shape_rate_to_constrained(in);
+    shape = temp(0);
+    rate = temp(1);
+  }
+
+  void set_from_proto(const ProtoState &state_, bool update_card) override {
+    if (update_card) {
+      card = state_.cardinality();
+    }
+    shape = state_.sr_state().shape();
+    rate = state_.sr_state().rate();
+  }
+
+  ProtoState get_as_proto() const override {
+    ProtoState state;
+    state.mutable_sr_state()->set_shape(shape);
+    state.mutable_sr_state()->set_rate(rate);
+    return state;
+  }
+
+  double log_det_jac() const override {
+    Eigen::VectorXd temp(2);
+    temp << shape, rate;
+    return shape_rate_log_det_jac(temp);
+  }
+};
+
+}  // namespace State
+
+#endif  // BAYESMIX_HIERARCHIES_LIKELIHOODS_STATES_SHAPE_RATE_STATE_H_

src/hierarchies/load_hierarchies.h

@@ -5,6 +5,7 @@
 #include <memory>
 
 #include "abstract_hierarchy.h"
+#include "betagg_hierarchy.h"
 #include "fa_hierarchy.h"
 #include "hierarchy_id.pb.h"
 #include "lapnig_hierarchy.h"
@@ -43,13 +44,17 @@ __attribute__((constructor)) static void load_hierarchies() {
   Builder<AbstractHierarchy> LapNIGbuilder = []() {
     return std::make_shared<LapNIGHierarchy>();
   };
+  Builder<AbstractHierarchy> BetaGGbuilder = []() {
+    return std::make_shared<BetaGGHierarchy>();
+  };
 
   factory.add_builder(NNIGHierarchy().get_id(), NNIGbuilder);
   factory.add_builder(NNxIGHierarchy().get_id(), NNxIGbuilder);
   factory.add_builder(NNWHierarchy().get_id(), NNWbuilder);
   factory.add_builder(LinRegUniHierarchy().get_id(), LinRegUnibuilder);
   factory.add_builder(FAHierarchy().get_id(), FAbuilder);
   factory.add_builder(LapNIGHierarchy().get_id(), LapNIGbuilder);
+  factory.add_builder(BetaGGHierarchy().get_id(), BetaGGbuilder);
 }
 
 #endif  // BAYESMIX_HIERARCHIES_LOAD_HIERARCHIES_H_

src/hierarchies/priors/CMakeLists.txt

@@ -13,4 +13,6 @@ target_sources(bayesmix PUBLIC
     mnig_prior_model.cc
     fa_prior_model.h
     fa_prior_model.cc
+    gamma_gamma_prior.h
+    gamma_gamma_prior.cc
 )