Feature/crispino nontransitive

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: BayesMallowsSMC2
 Type: Package
 Title: Nested Sequential Monte Carlo for the Bayesian Mallows Model
-Version: 0.3.0
+Version: 0.3.0.9000
 Authors@R: c(person("Oystein", "Sorensen",
                     email = "oystein.sorensen.1985@gmail.com",
                     role = c("aut", "cre"),

NEWS.md

@@ -1,3 +1,9 @@
+# BayesMallowsSMC2 (development version)
+
+## Major changes
+
+* Implemented the Bernoulli error model for non-transitive pairwise preferences (Crispino et al., 2019). The error probability `epsilon` is tracked and estimated. You can enable it using `error_model = "bernoulli"` in `set_smc_options()`, and configure the Beta prior via `kappa_1` and `kappa_2` in `set_hyperparameters()`.
+
 # BayesMallowsSMC2 version 0.3.0
 
 ## Major changes

R/set_hyperparameters.R

@@ -13,9 +13,13 @@
 #'   distribution for cluster probabilities. Only used when `n_clusters > 1`.
 #'   Defaults to 10.
 #' @param n_clusters Integer defining the number of clusters. Defaults to 1.
+#' @param kappa_1 First shape parameter of the Beta prior distribution for the
+#'   error probability epsilon. Defaults to 1.
+#' @param kappa_2 Second shape parameter of the Beta prior distribution for the
+#'   error probability epsilon. Defaults to 1.
 #'
 #' @return A list with components `n_items`, `alpha_shape`, `alpha_rate`,
-#'   `cluster_concentration`, and `n_clusters`.
+#'   `cluster_concentration`, `n_clusters`, `kappa_1`, and `kappa_2`.
 #' @export
 #'
 #' @examples
@@ -47,7 +51,7 @@
 #'
 set_hyperparameters <- function(
     n_items, alpha_shape = 1, alpha_rate = .5, cluster_concentration = 10,
-    n_clusters = 1) {
+    n_clusters = 1, kappa_1 = 1, kappa_2 = 1) {
   if(missing(n_items)) stop("n_items must be provided")
   as.list(environment())
 }

R/set_smc_options.R

@@ -50,6 +50,8 @@
 #' @param use_backward_simulation Logical specifying whether to use
 #'   Particle Gibbs with Backward Simulation (PG-BSi) during the rejuvenation
 #'   step. Defaults to `FALSE`.
+#' @param error_model Character string specifying the error model for pairwise
+#'   preferences. Options are `"none"` (default) or `"bernoulli"`.
 #'
 #' @details
 #' The SMC2 algorithm uses a nested particle filter structure:
@@ -129,6 +131,7 @@ set_smc_options <- function(
     max_rejuvenation_steps = 20,
     metric = "footrule", resampler = "multinomial",
     latent_rank_proposal = "uniform", verbose = FALSE,
-    trace = FALSE, trace_latent = FALSE, use_backward_simulation = FALSE) {
+    trace = FALSE, trace_latent = FALSE, use_backward_simulation = FALSE,
+    error_model = "none") {
   as.list(environment())
 }

src/options.cpp

@@ -13,4 +13,5 @@ Options::Options(const Rcpp::List& input_options) :
   verbose{input_options["verbose"]},
   trace{input_options["trace"]},
   trace_latent{input_options["trace_latent"]},
-  use_backward_simulation{input_options["use_backward_simulation"]}{}
+  use_backward_simulation{input_options["use_backward_simulation"]},
+  error_model ( input_options["error_model"] ){}

src/options.h

@@ -19,4 +19,5 @@ struct Options{
   const bool trace;
   const bool trace_latent;
   const bool use_backward_simulation;
+  const std::string error_model;
 };

src/parameter_tracer.cpp

@@ -27,6 +27,12 @@ void ParameterTracer::update_trace(const std::vector<Particle>& pvec, int t) {
     }
     tau_traces.push_back(tau);
 
+    vec epsilon(pvec.size());
+    for(size_t i{}; i < pvec.size(); i++) {
+      epsilon(i) = pvec[i].parameters.epsilon;
+    }
+    epsilon_traces.push_back(epsilon);
+
     vec log_importance_weights(pvec.size());
     for(size_t i{}; i < pvec.size(); i++) {
       log_importance_weights(i) = pvec[i].log_importance_weight;

src/parameter_tracer.h

@@ -10,6 +10,7 @@ struct ParameterTracer{
   std::vector<arma::mat> alpha_traces{};
   std::vector<arma::ucube> rho_traces{};
   std::vector<arma::mat> tau_traces{};
+  std::vector<arma::vec> epsilon_traces{};
   std::vector<arma::vec> log_importance_weights_traces{};
   std::vector<std::vector<arma::umat>> latent_rankings_traces{};
   void update_trace(const std::vector<Particle>& pvec, int t);

src/particle.cpp

@@ -7,14 +7,22 @@
 
 using namespace arma;
 
-StaticParameters::StaticParameters(const vec& alpha, const umat& rho, const vec& tau) :
-  alpha { alpha }, rho { rho }, tau { tau } {}
+StaticParameters::StaticParameters(const vec& alpha, const umat& rho, const vec& tau, double epsilon) :
+  alpha { alpha }, rho { rho }, tau { tau }, epsilon { epsilon } {}
 
-StaticParameters::StaticParameters(const Prior& prior) :
+StaticParameters::StaticParameters(const Prior& prior, const Options& options) :
   alpha { Rcpp::rgamma(prior.n_clusters, prior.alpha_shape, 1 / prior.alpha_rate) },
   rho { umat(prior.n_items, prior.n_clusters) },
   tau { normalise(Rcpp::as<vec>(Rcpp::rgamma(prior.n_clusters, prior.cluster_concentration, 1)), 1) }
   {
+    if (options.error_model == "bernoulli") {
+      double max_p = R::pbeta(0.5, prior.kappa_1, prior.kappa_2, 1, 0);
+      double u = R::runif(0, 1) * max_p;
+      epsilon = R::qbeta(u, prior.kappa_1, prior.kappa_2, 1, 0);
+      if (epsilon == 0.0) epsilon = 1e-6;
+    } else {
+      epsilon = 0.0;
+    }
     rho.each_col([&prior](uvec& a){
       a = Rcpp::as<uvec>(Rcpp::sample(prior.n_items, prior.n_items, false));
       });
@@ -57,7 +65,7 @@ void Particle::run_particle_filter(
   unsigned int pf_index{};
   for(auto& pf : particle_filters) {
     auto proposal = sample_latent_rankings(
-      data, t, prior, options.latent_rank_proposal, parameters, pfun, distfun);
+      data, t, prior, options.latent_rank_proposal, options.error_model, parameters, pfun, distfun);
 
     if(conditional && pf_index == 0) {
       if (options.use_backward_simulation) {
@@ -72,13 +80,41 @@ void Particle::run_particle_filter(
 
     double log_prob{};
 
-    for(size_t i{}; i < proposal.proposal.n_cols; i++) {
-      vec log_cluster_contribution(prior.n_clusters);
-      for(size_t c{}; c < prior.n_clusters; c++) {
-        log_cluster_contribution(c) = log(parameters.tau(c)) - this->logz(c) -
-          parameters.alpha(c) * distfun->d(proposal.proposal.col(i), parameters.rho.col(c));
+    if (options.error_model == "bernoulli" && dynamic_cast<PairwisePreferences*>(data.get())) {
+      PairwisePreferences* pp = dynamic_cast<PairwisePreferences*>(data.get());
+      pairwise_tp new_data = pp->timeseries[t];
+      size_t user_idx = 0;
+      for (auto ndit = new_data.begin(); ndit != new_data.end(); ++ndit) {
+        vec log_cluster_contribution(prior.n_clusters);
+        for(size_t c{}; c < prior.n_clusters; c++) {
+          log_cluster_contribution(c) = log(parameters.tau(c)) - this->logz(c) -
+            parameters.alpha(c) * distfun->d(proposal.proposal.col(user_idx), parameters.rho.col(c));
+        }
+        log_prob += log_sum_exp(log_cluster_contribution);
+
+        int a = 0;
+        int p_n = ndit->second.size();
+        for (auto pair : ndit->second) {
+          unsigned int item_A = pair.first - 1;
+          unsigned int item_B = pair.second - 1;
+          unsigned int rank_A = proposal.proposal(item_A, user_idx);
+          unsigned int rank_B = proposal.proposal(item_B, user_idx);
+          if (rank_A > rank_B) {
+            a++;
+          }
+        }
+        log_prob += a * log(parameters.epsilon / (1.0 - parameters.epsilon)) + p_n * log(1.0 - parameters.epsilon);
+        user_idx++;
+      }
+    } else {
+      for(size_t i{}; i < proposal.proposal.n_cols; i++) {
+        vec log_cluster_contribution(prior.n_clusters);
+        for(size_t c{}; c < prior.n_clusters; c++) {
+          log_cluster_contribution(c) = log(parameters.tau(c)) - this->logz(c) -
+            parameters.alpha(c) * distfun->d(proposal.proposal.col(i), parameters.rho.col(c));
+        }
+        log_prob += log_sum_exp(log_cluster_contribution);
       }
-      log_prob += log_sum_exp(log_cluster_contribution);
     }
 
     pf.cluster_probabilities = join_horiz(
@@ -120,7 +156,7 @@ std::vector<Particle> create_particle_vector(const Options& options, const Prior
   result.reserve(options.n_particles);
 
   for(size_t i{}; i < options.n_particles; i++) {
-    result.push_back(Particle{options, StaticParameters(prior), pfun});
+    result.push_back(Particle{options, StaticParameters(prior, options), pfun});
   }
 
   return result;

src/particle.h

@@ -10,11 +10,12 @@
 
 struct StaticParameters{
   StaticParameters() {}
-  StaticParameters(const arma::vec& alpha, const arma::umat& rho, const arma::vec& tau);
-  StaticParameters(const Prior& prior);
+  StaticParameters(const arma::vec& alpha, const arma::umat& rho, const arma::vec& tau, double epsilon);
+  StaticParameters(const Prior& prior, const Options& options);
   arma::vec alpha;
   arma::umat rho;
   arma::vec tau;
+  double epsilon;
 };
 
 struct ParticleFilter{

src/prior.cpp

@@ -5,4 +5,6 @@ Prior::Prior(const Rcpp::List& input_prior) :
   alpha_rate { input_prior["alpha_rate"] },
   cluster_concentration { input_prior["cluster_concentration"] },
   n_clusters { input_prior["n_clusters"] },
-  n_items { input_prior["n_items"]} {}
+  n_items { input_prior["n_items"]},
+  kappa_1 { input_prior.containsElementNamed("kappa_1") ? Rcpp::as<double>(input_prior["kappa_1"]) : 1.0 },
+  kappa_2 { input_prior.containsElementNamed("kappa_2") ? Rcpp::as<double>(input_prior["kappa_2"]) : 1.0 } {}

src/prior.h

@@ -9,4 +9,6 @@ struct Prior{
   int cluster_concentration;
   int n_clusters;
   int n_items;
+  double kappa_1;
+  double kappa_2;
 };

src/rejuvenate.cpp

@@ -68,7 +68,7 @@ bool Particle::rejuvenate(
     rho_proposal.col(cluster) = leap_and_shift(parameters.rho.col(cluster), cluster, prior);
   }
 
-  Particle proposal_particle(options, StaticParameters{alpha_proposal, rho_proposal, parameters.tau}, pfun);
+  Particle proposal_particle(options, StaticParameters{alpha_proposal, rho_proposal, parameters.tau, parameters.epsilon}, pfun);
 
   double log_ratio{};
   vec additional_terms = prior.alpha_shape * (log(alpha_proposal) - log(parameters.alpha)) -
@@ -86,7 +86,7 @@ bool Particle::rejuvenate(
 
   bool accepted{};
   if(log_ratio > log(R::runif(0, 1))) {
-    this->parameters = StaticParameters{alpha_proposal, rho_proposal, parameters.tau};
+    this->parameters = StaticParameters{alpha_proposal, rho_proposal, parameters.tau, parameters.epsilon};
     this->conditioned_particle_filter = proposed_particle_filter;
     this->log_incremental_likelihood = proposal_particle.log_incremental_likelihood;
     this->log_normalized_particle_filter_weights = proposal_particle.log_normalized_particle_filter_weights;
@@ -97,14 +97,17 @@ bool Particle::rejuvenate(
     accepted = false;
   }
 
-  if(prior.n_clusters > 1) {
-    uvec cluster_assignments = particle_filters[conditioned_particle_filter].cluster_assignments;
-    uvec cluster_frequencies = hist(cluster_assignments, regspace<uvec>(0, prior.n_clusters - 1));
+  if(prior.n_clusters > 1 || options.error_model == "bernoulli") {
+    if(prior.n_clusters > 1) {
+      uvec cluster_assignments = particle_filters[conditioned_particle_filter].cluster_assignments;
+      uvec cluster_frequencies = hist(cluster_assignments, regspace<uvec>(0, prior.n_clusters - 1));
 
-    for(size_t cluster{}; cluster < prior.n_clusters; cluster++) {
-      parameters.tau(cluster) = R::rgamma(cluster_frequencies(cluster) + prior.cluster_concentration, 1.0);
+      for(size_t cluster{}; cluster < prior.n_clusters; cluster++) {
+        parameters.tau(cluster) = R::rgamma(cluster_frequencies(cluster) + prior.cluster_concentration, 1.0);
+      }
+      parameters.tau = normalise(parameters.tau, 1);
     }
-    parameters.tau = normalise(parameters.tau, 1);
+
     Particle gibbs_particle(options, this->parameters, pfun);
     gibbs_particle.conditioned_particle_filter = 0;
 
@@ -122,15 +125,50 @@ bool Particle::rejuvenate(
         int b_t = Rcpp::sample(probs.size(), 1, false, probs, false)[0];
 
         bsi_latent_rankings.col(t) = this->particle_filters[b_t].latent_rankings.col(t);
-        bsi_cluster_assignments(t) = this->particle_filters[b_t].cluster_assignments(t);
+        if (prior.n_clusters > 1) {
+          bsi_cluster_assignments(t) = this->particle_filters[b_t].cluster_assignments(t);
+        }
       }
       gibbs_particle.reference_latent_rankings = bsi_latent_rankings;
-      gibbs_particle.reference_cluster_assignments = bsi_cluster_assignments;
+      if (prior.n_clusters > 1) {
+        gibbs_particle.reference_cluster_assignments = bsi_cluster_assignments;
+      }
     } else {
       gibbs_particle.reference_latent_rankings = this->particle_filters[this->conditioned_particle_filter].latent_rankings;
-      gibbs_particle.reference_cluster_assignments = this->particle_filters[this->conditioned_particle_filter].cluster_assignments;
+      if (prior.n_clusters > 1) {
+        gibbs_particle.reference_cluster_assignments = this->particle_filters[this->conditioned_particle_filter].cluster_assignments;
+      }
     }
 
+    if (options.error_model == "bernoulli" && dynamic_cast<PairwisePreferences*>(data.get())) {
+      int total_a = 0;
+      int total_p = 0;
+      PairwisePreferences* pp = dynamic_cast<PairwisePreferences*>(data.get());
+      for (size_t t{}; t < T + 1; t++) {
+        const pairwise_tp& new_data = pp->timeseries[t];
+        for (auto ndit = new_data.begin(); ndit != new_data.end(); ++ndit) {
+          total_p += ndit->second.size();
+          for (const auto& pair : ndit->second) {
+            unsigned int item_A = pair.first - 1;
+            unsigned int item_B = pair.second - 1;
+            unsigned int rank_A = gibbs_particle.reference_latent_rankings(item_A, t);
+            unsigned int rank_B = gibbs_particle.reference_latent_rankings(item_B, t);
+            if (rank_A > rank_B) {
+              total_a++;
+            }
+          }
+        }
+      }
+      int total_b = total_p - total_a;
+
+      double max_p = R::pbeta(0.5, prior.kappa_1 + total_a, prior.kappa_2 + total_b, 1, 0);
+      double u = R::runif(0, 1) * max_p;
+      double epsilon_prime = R::qbeta(u, prior.kappa_1 + total_a, prior.kappa_2 + total_b, 1, 0);
+      if (epsilon_prime == 0.0) epsilon_prime = 1e-6;
+      parameters.epsilon = epsilon_prime;
+      gibbs_particle.parameters.epsilon = epsilon_prime;
+    }
+
     if (options.use_backward_simulation) {
       gibbs_particle.particle_filters[0] = ParticleFilter{};
       gibbs_particle.particle_filters[0].cluster_probabilities = mat{};

src/run_smc.cpp

@@ -115,6 +115,7 @@ Rcpp::List run_smc(
   mat alpha(prior.n_clusters, particle_vector.size());
   ucube rho(prior.n_items, prior.n_clusters, particle_vector.size());
   mat tau(prior.n_clusters, particle_vector.size());
+  vec epsilon(particle_vector.size());
   cube cluster_probabilities;
   if(prior.n_clusters > 1) {
     cluster_probabilities = cube(particle_vector.size(), particle_vector[0].particle_filters[0].cluster_probabilities.n_cols, prior.n_clusters);
@@ -124,6 +125,7 @@ Rcpp::List run_smc(
     alpha.col(i) = particle_vector[i].parameters.alpha;
     rho.slice(i) = particle_vector[i].parameters.rho;
     tau.col(i) = particle_vector[i].parameters.tau;
+    epsilon(i) = particle_vector[i].parameters.epsilon;
 
     if(prior.n_clusters > 1) {
       cluster_probabilities.row(i) = particle_vector[i].particle_filters[particle_vector[i].conditioned_particle_filter].cluster_probabilities.t();
@@ -134,6 +136,7 @@ Rcpp::List run_smc(
     Rcpp::Named("alpha") = alpha,
     Rcpp::Named("rho") = rho,
     Rcpp::Named("tau") = tau,
+    Rcpp::Named("epsilon") = epsilon,
     Rcpp::Named("cluster_probabilities") = cluster_probabilities,
     Rcpp::Named("ESS") = ESS,
     Rcpp::Named("resampling") = resampling,
@@ -143,6 +146,7 @@ Rcpp::List run_smc(
     Rcpp::Named("alpha_traces") = tracer.alpha_traces,
     Rcpp::Named("rho_traces") = tracer.rho_traces,
     Rcpp::Named("tau_traces") = tracer.tau_traces,
+    Rcpp::Named("epsilon_traces") = tracer.epsilon_traces,
     Rcpp::Named("log_importance_weights_traces") = tracer.log_importance_weights_traces,
     Rcpp::Named("latent_rankings_traces") = tracer.latent_rankings_traces
   );