From e069f9ad0c7aa54eb38bb6c68f5d624193c0bef7 Mon Sep 17 00:00:00 2001
From: Chris Dyer <cdyer@cs.cmu.edu>
Date: Tue, 3 Jan 2012 16:59:11 -0500
Subject: multi-floor chinese restaurant described by wood&teh (2009)
---
utils/mfcr.h | 354 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 354 insertions(+)
create mode 100644 utils/mfcr.h
(limited to 'utils/mfcr.h')
diff --git a/utils/mfcr.h b/utils/mfcr.h
new file mode 100644
index 00000000..3eb133fc
--- /dev/null
+++ b/utils/mfcr.h
@@ -0,0 +1,354 @@
+#ifndef _MFCR_H_
+#define _MFCR_H_
+
+#include <algorithm>
+#include <numeric>
+#include <cassert>
+#include <cmath>
+#include <list>
+#include <iostream>
+#include <vector>
+#include <tr1/unordered_map>
+#include <boost/functional/hash.hpp>
+#include "sampler.h"
+#include "slice_sampler.h"
+
+struct TableCount {
+ TableCount() : count(), floor() {}
+ TableCount(int c, int f) : count(c), floor(f) {
+ assert(f >= 0);
+ }
+ int count; // count or delta (may be 0, <0, or >0)
+ unsigned char floor; // from which floor?
+};
+
+std::ostream& operator<<(std::ostream& o, const TableCount& tc) {
+ return o << "[c=" << tc.count << " floor=" << static_cast<unsigned int>(tc.floor) << ']';
+}
+
+// Multi-Floor Chinese Restaurant as proposed by Wood & Teh (AISTATS, 2009) to simulate
+// graphical Pitman-Yor processes.
+// http://jmlr.csail.mit.edu/proceedings/papers/v5/wood09a/wood09a.pdf
+//
+// Implementation is based on Blunsom, Cohn, Goldwater, & Johnson (ACL 2009) and code
+// referenced therein.
+// http://www.aclweb.org/anthology/P/P09/P09-2085.pdf
+//
+template <typename Dish, typename DishHash = boost::hash<Dish> >
+class MFCR {
+ public:
+
+ MFCR(unsigned num_floors, double d, double alpha) :
+ num_floors_(num_floors),
+ num_tables_(),
+ num_customers_(),
+ d_(d),
+ alpha_(alpha),
+ d_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
+ d_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
+ alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+ alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+
+ MFCR(unsigned num_floors, double d_alpha, double d_beta, double alpha_shape, double alpha_rate, double d = 0.9, double alpha = 10.0) :
+ num_floors_(num_floors),
+ num_tables_(),
+ num_customers_(),
+ d_(d),
+ alpha_(alpha),
+ d_prior_alpha_(d_alpha),
+ d_prior_beta_(d_beta),
+ alpha_prior_shape_(alpha_shape),
+ alpha_prior_rate_(alpha_rate) {}
+
+ double d() const { return d_; }
+ double alpha() const { return alpha_; }
+
+ bool has_d_prior() const {
+ return !std::isnan(d_prior_alpha_);
+ }
+
+ bool has_alpha_prior() const {
+ return !std::isnan(alpha_prior_shape_);
+ }
+
+ void clear() {
+ num_tables_ = 0;
+ num_customers_ = 0;
+ dish_locs_.clear();
+ }
+
+ unsigned num_tables() const {
+ return num_tables_;
+ }
+
+ unsigned num_tables(const Dish& dish) const {
+ const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
+ if (it == dish_locs_.end()) return 0;
+ return it->second.table_counts_.size();
+ }
+
+ // this is not terribly efficient but it should not typically be necessary to execute this query
+ unsigned num_tables(const Dish& dish, const unsigned floor) const {
+ const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
+ if (it == dish_locs_.end()) return 0;
+ unsigned c = 0;
+ for (typename std::list<TableCount>::const_iterator i = it->second.table_counts_.begin();
+ i != it->second.table_counts_.end(); ++i) {
+ if (i->floor == floor) ++c;
+ }
+ return c;
+ }
+
+ unsigned num_customers() const {
+ return num_customers_;
+ }
+
+ unsigned num_customers(const Dish& dish) const {
+ const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
+ if (it == dish_locs_.end()) return 0;
+ return it->total_dish_count_;
+ }
+
+ // returns (delta, floor) indicating whether a new table (delta) was opened and on which floor
+ TableCount increment(const Dish& dish, const std::vector<double>& p0s, const std::vector<double>& lambdas, MT19937* rng) {
+ assert(p0s.size() == num_floors_);
+ assert(lambdas.size() == num_floors_);
+
+ DishLocations& loc = dish_locs_[dish];
+ // marg_p0 = marginal probability of opening a new table on any floor with label dish
+ const double marg_p0 = std::inner_product(p0s.begin(), p0s.end(), lambdas.begin(), 0.0);
+ assert(marg_p0 <= 1.0);
+ int floor = -1;
+ bool share_table = false;
+ if (loc.total_dish_count_) {
+ const double p_empty = (alpha_ + num_tables_ * d_) * marg_p0;
+ const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * d_);
+ share_table = rng->SelectSample(p_empty, p_share);
+ }
+ if (share_table) {
+ double r = rng->next() * (loc.total_dish_count_ - loc.table_counts_.size() * d_);
+ for (typename std::list<TableCount>::iterator ti = loc.table_counts_.begin();
+ ti != loc.table_counts_.end(); ++ti) {
+ r -= ti->count - d_;
+ if (r <= 0.0) {
+ ++ti->count;
+ floor = ti->floor;
+ break;
+ }
+ }
+ if (r > 0.0) {
+ std::cerr << "Serious error: r=" << r << std::endl;
+ Print(&std::cerr);
+ assert(r <= 0.0);
+ }
+ } else { // sit at currently empty table -- must sample what floor
+ double r = rng->next() * marg_p0;
+ for (unsigned i = 0; i < p0s.size(); ++i) {
+ r -= p0s[i] * lambdas[i];
+ if (r <= 0.0) {
+ floor = i;
+ break;
+ }
+ }
+ assert(floor >= 0);
+ loc.table_counts_.push_back(TableCount(1, floor));
+ ++num_tables_;
+ }
+ ++loc.total_dish_count_;
+ ++num_customers_;
+ return (share_table ? TableCount(0, floor) : TableCount(1, floor));
+ }
+
+ // returns first = -1 or 0, indicating whether a table was closed, and on what floor (second)
+ TableCount decrement(const Dish& dish, MT19937* rng) {
+ DishLocations& loc = dish_locs_[dish];
+ assert(loc.total_dish_count_);
+ int floor = -1;
+ int delta = 0;
+ if (loc.total_dish_count_ == 1) {
+ floor = loc.table_counts_.front().floor;
+ dish_locs_.erase(dish);
+ --num_tables_;
+ --num_customers_;
+ delta = -1;
+ } else {
+ // sample customer to remove UNIFORMLY. that is, do NOT use the d
+ // here. if you do, it will introduce (unwanted) bias!
+ double r = rng->next() * loc.total_dish_count_;
+ --loc.total_dish_count_;
+ --num_customers_;
+ for (typename std::list<TableCount>::iterator ti = loc.table_counts_.begin();
+ ti != loc.table_counts_.end(); ++ti) {
+ r -= ti->count;
+ if (r <= 0.0) {
+ floor = ti->floor;
+ if ((--ti->count) == 0) {
+ --num_tables_;
+ delta = -1;
+ loc.table_counts_.erase(ti);
+ }
+ break;
+ }
+ }
+ if (r > 0.0) {
+ std::cerr << "Serious error: r=" << r << std::endl;
+ Print(&std::cerr);
+ assert(r <= 0.0);
+ }
+ }
+ return TableCount(delta, floor);
+ }
+
+ double prob(const Dish& dish, const std::vector<double>& p0s, const std::vector<double>& lambdas) const {
+ assert(p0s.size() == num_floors_);
+ assert(lambdas.size() == num_floors_);
+ const double marg_p0 = std::inner_product(p0s.begin(), p0s.end(), lambdas.begin(), 0.0);
+ assert(marg_p0 <= 1.0);
+ const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
+ const double r = num_tables_ * d_ + alpha_;
+ if (it == dish_locs_.end()) {
+ return r * marg_p0 / (num_customers_ + alpha_);
+ } else {
+ return (it->second.total_dish_count_ - d_ * it->second.table_counts_.size() + r * marg_p0) /
+ (num_customers_ + alpha_);
+ }
+ }
+
+ double log_crp_prob() const {
+ return log_crp_prob(d_, alpha_);
+ }
+
+ static double log_beta_density(const double& x, const double& alpha, const double& beta) {
+ assert(x > 0.0);
+ assert(x < 1.0);
+ assert(alpha > 0.0);
+ assert(beta > 0.0);
+ const double lp = (alpha-1)*log(x)+(beta-1)*log(1-x)+lgamma(alpha+beta)-lgamma(alpha)-lgamma(beta);
+ return lp;
+ }
+
+ static double log_gamma_density(const double& x, const double& shape, const double& rate) {
+ assert(x >= 0.0);
+ assert(shape > 0.0);
+ assert(rate > 0.0);
+ const double lp = (shape-1)*log(x) - shape*log(rate) - x/rate - lgamma(shape);
+ return lp;
+ }
+
+ // taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
+ // does not include draws from G_w's
+ double log_crp_prob(const double& d, const double& alpha) const {
+ double lp = 0.0;
+ if (has_d_prior())
+ lp = log_beta_density(d, d_prior_alpha_, d_prior_beta_);
+ if (has_alpha_prior())
+ lp += log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
+ assert(lp <= 0.0);
+ if (num_customers_) {
+ if (d > 0.0) {
+ const double r = lgamma(1.0 - d);
+ lp += lgamma(alpha) - lgamma(alpha + num_customers_)
+ + num_tables_ * log(d) + lgamma(alpha / d + num_tables_)
+ - lgamma(alpha / d);
+ assert(std::isfinite(lp));
+ for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
+ it != dish_locs_.end(); ++it) {
+ const DishLocations& cur = it->second;
+ for (std::list<TableCount>::const_iterator ti = cur.table_counts_.begin(); ti != cur.table_counts_.end(); ++ti) {
+ lp += lgamma(ti->count - d) - r;
+ }
+ }
+ } else {
+ assert(!"not implemented yet");
+ }
+ }
+ assert(std::isfinite(lp));
+ return lp;
+ }
+
+ void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
+ assert(has_d_prior() || has_alpha_prior());
+ DiscountResampler dr(*this);
+ ConcentrationResampler cr(*this);
+ for (int iter = 0; iter < nloop; ++iter) {
+ if (has_alpha_prior()) {
+ alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
+ std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+ }
+ if (has_d_prior()) {
+ d_ = slice_sampler1d(dr, d_, *rng, std::numeric_limits<double>::min(),
+ 1.0, 0.0, niterations, 100*niterations);
+ }
+ }
+ alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
+ std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+ }
+
+ struct DiscountResampler {
+ DiscountResampler(const MFCR& crp) : crp_(crp) {}
+ const MFCR& crp_;
+ double operator()(const double& proposed_d) const {
+ return crp_.log_crp_prob(proposed_d, crp_.alpha_);
+ }
+ };
+
+ struct ConcentrationResampler {
+ ConcentrationResampler(const MFCR& crp) : crp_(crp) {}
+ const MFCR& crp_;
+ double operator()(const double& proposed_alpha) const {
+ return crp_.log_crp_prob(crp_.d_, proposed_alpha);
+ }
+ };
+
+ struct DishLocations {
+ DishLocations() : total_dish_count_() {}
+ unsigned total_dish_count_; // customers at all tables with this dish
+ std::list<TableCount> table_counts_; // list<> gives O(1) deletion and insertion, which we want
+ // .size() is the number of tables for this dish
+ };
+
+ void Print(std::ostream* out) const {
+ (*out) << "MFCR(d=" << d_ << ",alpha=" << alpha_ << ") customers=" << num_customers_ << std::endl;
+ for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
+ it != dish_locs_.end(); ++it) {
+ (*out) << it->first << " (" << it->second.total_dish_count_ << " on " << it->second.table_counts_.size() << " tables): ";
+ for (typename std::list<TableCount>::const_iterator i = it->second.table_counts_.begin();
+ i != it->second.table_counts_.end(); ++i) {
+ (*out) << " " << *i;
+ }
+ (*out) << std::endl;
+ }
+ }
+
+ typedef typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator const_iterator;
+ const_iterator begin() const {
+ return dish_locs_.begin();
+ }
+ const_iterator end() const {
+ return dish_locs_.end();
+ }
+
+ unsigned num_floors_;
+ unsigned num_tables_;
+ unsigned num_customers_;
+ std::tr1::unordered_map<Dish, DishLocations, DishHash> dish_locs_;
+
+ double d_;
+ double alpha_;
+
+ // optional beta prior on d_ (NaN if no prior)
+ double d_prior_alpha_;
+ double d_prior_beta_;
+
+ // optional gamma prior on alpha_ (NaN if no prior)
+ double alpha_prior_shape_;
+ double alpha_prior_rate_;
+};
+
+template <typename T,typename H>
+std::ostream& operator<<(std::ostream& o, const MFCR<T,H>& c) {
+ c.Print(&o);
+ return o;
+}
+
+#endif
--
cgit v1.2.3
From 400d60b20e9e480b0eff9843404a4cb9f8bd02cc Mon Sep 17 00:00:00 2001
From: Chris Dyer <cdyer@cs.cmu.edu>
Date: Wed, 8 Feb 2012 16:22:55 -0500
Subject: move widely duplicated math functions into m.h header
---
.gitignore | 1 +
gi/pf/base_distributions.cc | 22 +++++------
gi/pf/base_distributions.h | 21 +---------
gi/pf/conditional_pseg.h | 3 +-
gi/pf/pfdist.cc | 6 +--
gi/pf/pfnaive.cc | 4 +-
phrasinator/gibbs_train_plm.cc | 8 +---
utils/Makefile.am | 5 ++-
utils/m.h | 89 ++++++++++++++++++++++++++++++++++++++++++
utils/m_test.cc | 75 +++++++++++++++++++++++++++++++++++
utils/mfcr.h | 22 ++---------
11 files changed, 194 insertions(+), 62 deletions(-)
create mode 100644 utils/m.h
create mode 100644 utils/m_test.cc
(limited to 'utils/mfcr.h')
diff --git a/.gitignore b/.gitignore
index ab8bf2c7..4f75d153 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
mira/kbest_mira
+utils/m_test
sa-extract/calignment.c
sa-extract/calignment.so
sa-extract/cdat.c
diff --git a/gi/pf/base_distributions.cc b/gi/pf/base_distributions.cc
index d362fd76..d9761005 100644
--- a/gi/pf/base_distributions.cc
+++ b/gi/pf/base_distributions.cc
@@ -59,7 +59,7 @@ prob_t PhraseConditionalUninformativeUnigramBase::p0(const vector<WordID>& vsrc,
const int flen = vsrc.size() - start_src;
const int elen = vtrg.size() - start_trg;
prob_t p;
- p.logeq(log_poisson(elen, flen + 0.01)); // elen | flen ~Pois(flen + 0.01)
+ p.logeq(Md::log_poisson(elen, flen + 0.01)); // elen | flen ~Pois(flen + 0.01)
//p.logeq(log_poisson(elen, 1)); // elen | flen ~Pois(flen + 0.01)
for (int i = 0; i < elen; ++i)
p *= u(vtrg[i + start_trg]); // draw e_i ~Uniform
@@ -73,7 +73,7 @@ prob_t PhraseConditionalUninformativeBase::p0(const vector<WordID>& vsrc,
const int elen = vtrg.size() - start_trg;
prob_t p;
//p.logeq(log_poisson(elen, flen + 0.01)); // elen | flen ~Pois(flen + 0.01)
- p.logeq(log_poisson(elen, 1)); // elen | flen ~Pois(flen + 0.01)
+ p.logeq(Md::log_poisson(elen, 1)); // elen | flen ~Pois(flen + 0.01)
for (int i = 0; i < elen; ++i)
p *= kUNIFORM_TARGET; // draw e_i ~Uniform
return p;
@@ -113,7 +113,7 @@ prob_t PhraseConditionalBase::p0(const vector<WordID>& vsrc,
const int elen = vtrg.size() - start_trg;
prob_t uniform_src_alignment; uniform_src_alignment.logeq(-log(flen + 1));
prob_t p;
- p.logeq(log_poisson(elen, flen + 0.01)); // elen | flen ~Pois(flen + 0.01)
+ p.logeq(Md::log_poisson(elen, flen + 0.01)); // elen | flen ~Pois(flen + 0.01)
for (int i = 0; i < elen; ++i) { // for each position i in e-RHS
const WordID trg = vtrg[i + start_trg];
prob_t tp = prob_t::Zero();
@@ -139,9 +139,9 @@ prob_t PhraseJointBase::p0(const vector<WordID>& vsrc,
const int elen = vtrg.size() - start_trg;
prob_t uniform_src_alignment; uniform_src_alignment.logeq(-log(flen + 1));
prob_t p;
- p.logeq(log_poisson(flen, 1.0)); // flen ~Pois(1)
+ p.logeq(Md::log_poisson(flen, 1.0)); // flen ~Pois(1)
// elen | flen ~Pois(flen + 0.01)
- prob_t ptrglen; ptrglen.logeq(log_poisson(elen, flen + 0.01));
+ prob_t ptrglen; ptrglen.logeq(Md::log_poisson(elen, flen + 0.01));
p *= ptrglen;
p *= kUNIFORM_SOURCE.pow(flen); // each f in F ~Uniform
for (int i = 0; i < elen; ++i) { // for each position i in E
@@ -171,9 +171,9 @@ prob_t PhraseJointBase_BiDir::p0(const vector<WordID>& vsrc,
prob_t uniform_trg_alignment; uniform_trg_alignment.logeq(-log(elen + 1));
prob_t p1;
- p1.logeq(log_poisson(flen, 1.0)); // flen ~Pois(1)
+ p1.logeq(Md::log_poisson(flen, 1.0)); // flen ~Pois(1)
// elen | flen ~Pois(flen + 0.01)
- prob_t ptrglen; ptrglen.logeq(log_poisson(elen, flen + 0.01));
+ prob_t ptrglen; ptrglen.logeq(Md::log_poisson(elen, flen + 0.01));
p1 *= ptrglen;
p1 *= kUNIFORM_SOURCE.pow(flen); // each f in F ~Uniform
for (int i = 0; i < elen; ++i) { // for each position i in E
@@ -193,9 +193,9 @@ prob_t PhraseJointBase_BiDir::p0(const vector<WordID>& vsrc,
}
prob_t p2;
- p2.logeq(log_poisson(elen, 1.0)); // elen ~Pois(1)
+ p2.logeq(Md::log_poisson(elen, 1.0)); // elen ~Pois(1)
// flen | elen ~Pois(flen + 0.01)
- prob_t psrclen; psrclen.logeq(log_poisson(flen, elen + 0.01));
+ prob_t psrclen; psrclen.logeq(Md::log_poisson(flen, elen + 0.01));
p2 *= psrclen;
p2 *= kUNIFORM_TARGET.pow(elen); // each f in F ~Uniform
for (int i = 0; i < flen; ++i) { // for each position i in E
@@ -227,9 +227,9 @@ JumpBase::JumpBase() : p(200) {
for (int j = min_jump; j <= max_jump; ++j) {
prob_t& cp = cpd[j];
if (j < 0)
- cp.logeq(log_poisson(1.5-j, 1));
+ cp.logeq(Md::log_poisson(1.5-j, 1));
else if (j > 0)
- cp.logeq(log_poisson(j, 1));
+ cp.logeq(Md::log_poisson(j, 1));
cp.poweq(0.2);
z += cp;
}
diff --git a/gi/pf/base_distributions.h b/gi/pf/base_distributions.h
index a23ac32b..0d597c5c 100644
--- a/gi/pf/base_distributions.h
+++ b/gi/pf/base_distributions.h
@@ -13,24 +13,7 @@
#include "prob.h"
#include "tdict.h"
#include "sampler.h"
-
-inline double log_poisson(unsigned x, const double& lambda) {
- assert(lambda > 0.0);
- return log(lambda) * x - lgamma(x + 1) - lambda;
-}
-
-inline double log_binom_coeff(unsigned n, unsigned k) {
- assert(n >= k);
- if (n == k) return 0.0;
- return lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1);
-}
-
-// http://en.wikipedia.org/wiki/Negative_binomial_distribution
-inline double log_negative_binom(unsigned x, unsigned r, double p) {
- assert(p > 0.0);
- assert(p < 1.0);
- return log_binom_coeff(x + r - 1, x) + r * log(1 - p) + x * log(p);
-}
+#include "m.h"
inline std::ostream& operator<<(std::ostream& os, const std::vector<WordID>& p) {
os << '[';
@@ -68,7 +51,7 @@ struct Model1 {
struct PoissonUniformUninformativeBase {
explicit PoissonUniformUninformativeBase(const unsigned ves) : kUNIFORM(1.0 / ves) {}
prob_t operator()(const TRule& r) const {
- prob_t p; p.logeq(log_poisson(r.e_.size(), 1.0));
+ prob_t p; p.logeq(Md::log_poisson(r.e_.size(), 1.0));
prob_t q = kUNIFORM; q.poweq(r.e_.size());
p *= q;
return p;
diff --git a/gi/pf/conditional_pseg.h b/gi/pf/conditional_pseg.h
index 0aa5e8e0..2e9e38fc 100644
--- a/gi/pf/conditional_pseg.h
+++ b/gi/pf/conditional_pseg.h
@@ -6,6 +6,7 @@
#include <boost/functional/hash.hpp>
#include <iostream>
+#include "m.h"
#include "prob.h"
#include "ccrp_nt.h"
#include "mfcr.h"
@@ -210,7 +211,7 @@ struct ConditionalParallelSegementationModel {
prob_t AlignProbability(unsigned span) const {
prob_t p;
- p.logeq(aligns.logprob(span, log_poisson(span, 1.0)));
+ p.logeq(aligns.logprob(span, Md::log_poisson(span, 1.0)));
return p;
}
diff --git a/gi/pf/pfdist.cc b/gi/pf/pfdist.cc
index ef08a165..3d578db2 100644
--- a/gi/pf/pfdist.cc
+++ b/gi/pf/pfdist.cc
@@ -315,7 +315,7 @@ struct BackwardEstimate {
for (int i = 0; i < src_cov.size(); ++i)
if (!src_cov[i]) r.push_back(src_[i]);
const prob_t uniform_alignment(1.0 / r.size());
- e.logeq(log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
+ e.logeq(Md::log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
for (unsigned j = trg_cov; j < trg_.size(); ++j) {
prob_t p;
for (unsigned i = 0; i < r.size(); ++i)
@@ -352,7 +352,7 @@ struct BackwardEstimateSym {
if (!src_cov[i]) r.push_back(src_[i]);
r.push_back(0); // NULL word
const prob_t uniform_alignment(1.0 / r.size());
- e.logeq(log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
+ e.logeq(Md::log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
for (unsigned j = trg_cov; j < trg_.size(); ++j) {
prob_t p;
for (unsigned i = 0; i < r.size(); ++i)
@@ -367,7 +367,7 @@ struct BackwardEstimateSym {
r.pop_back();
const prob_t inv_uniform(1.0 / (trg_.size() - trg_cov + 1.0));
prob_t inv;
- inv.logeq(log_poisson(r.size(), trg_.size() - trg_cov));
+ inv.logeq(Md::log_poisson(r.size(), trg_.size() - trg_cov));
for (unsigned i = 0; i < r.size(); ++i) {
prob_t p;
for (unsigned j = trg_cov - 1; j < trg_.size(); ++j)
diff --git a/gi/pf/pfnaive.cc b/gi/pf/pfnaive.cc
index acba9d22..e1a53f5c 100644
--- a/gi/pf/pfnaive.cc
+++ b/gi/pf/pfnaive.cc
@@ -77,7 +77,7 @@ struct BackwardEstimateSym {
r.push_back(src_[i]);
r.push_back(0); // NULL word
const prob_t uniform_alignment(1.0 / r.size());
- e.logeq(log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
+ e.logeq(Md::log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
for (unsigned j = trg_cov; j < trg_.size(); ++j) {
prob_t p;
for (unsigned i = 0; i < r.size(); ++i)
@@ -92,7 +92,7 @@ struct BackwardEstimateSym {
r.pop_back();
const prob_t inv_uniform(1.0 / (trg_.size() - trg_cov + 1.0));
prob_t inv;
- inv.logeq(log_poisson(r.size(), trg_.size() - trg_cov));
+ inv.logeq(Md::log_poisson(r.size(), trg_.size() - trg_cov));
for (unsigned i = 0; i < r.size(); ++i) {
prob_t p;
for (unsigned j = trg_cov - 1; j < trg_.size(); ++j)
diff --git a/phrasinator/gibbs_train_plm.cc b/phrasinator/gibbs_train_plm.cc
index 29b3d7ea..66b46011 100644
--- a/phrasinator/gibbs_train_plm.cc
+++ b/phrasinator/gibbs_train_plm.cc
@@ -8,6 +8,7 @@
#include "dict.h"
#include "sampler.h"
#include "ccrp.h"
+#include "m.h"
using namespace std;
using namespace std::tr1;
@@ -95,11 +96,6 @@ void ReadCorpus(const string& filename, vector<vector<int> >* c, set<int>* vocab
if (in != &cin) delete in;
}
-double log_poisson(unsigned x, const double& lambda) {
- assert(lambda > 0.0);
- return log(lambda) * x - lgamma(x + 1) - lambda;
-}
-
struct UniphraseLM {
UniphraseLM(const vector<vector<int> >& corpus,
const set<int>& vocab,
@@ -128,7 +124,7 @@ struct UniphraseLM {
double log_p0(const vector<int>& phrase) const {
double len_logprob;
if (use_poisson_)
- len_logprob = log_poisson(phrase.size(), 1.0);
+ len_logprob = Md::log_poisson(phrase.size(), 1.0);
else
len_logprob = log(1 - p_end_) * (phrase.size() -1) + log(p_end_);
return log(uniform_word_) * phrase.size() + len_logprob;
diff --git a/utils/Makefile.am b/utils/Makefile.am
index 3e559c75..a1ea8270 100644
--- a/utils/Makefile.am
+++ b/utils/Makefile.am
@@ -7,11 +7,12 @@ TESTS = ts phmt mfcr_test
if HAVE_GTEST
noinst_PROGRAMS += \
dict_test \
+ m_test \
weights_test \
logval_test \
small_vector_test
-TESTS += small_vector_test logval_test weights_test dict_test
+TESTS += small_vector_test logval_test weights_test dict_test m_test
endif
reconstruct_weights_SOURCES = reconstruct_weights.cc
@@ -38,6 +39,8 @@ endif
phmt_SOURCES = phmt.cc
ts_SOURCES = ts.cc
+m_test_SOURCES = m_test.cc
+m_test_LDADD = $(GTEST_LDFLAGS) $(GTEST_LIBS)
dict_test_SOURCES = dict_test.cc
dict_test_LDADD = $(GTEST_LDFLAGS) $(GTEST_LIBS)
mfcr_test_SOURCES = mfcr_test.cc
diff --git a/utils/m.h b/utils/m.h
new file mode 100644
index 00000000..b25248c2
--- /dev/null
+++ b/utils/m.h
@@ -0,0 +1,89 @@
+#ifndef _M_H_
+#define _M_H_
+
+#include <cassert>
+#include <cmath>
+
+template <typename F>
+struct M {
+ // support [0, 1, 2 ...)
+ static inline F log_poisson(unsigned x, const F& lambda) {
+ assert(lambda > 0.0);
+ return std::log(lambda) * x - lgamma(x + 1) - lambda;
+ }
+
+ // support [0, 1, 2 ...)
+ static inline F log_geometric(unsigned x, const F& p) {
+ assert(p > 0.0);
+ assert(p < 1.0);
+ return std::log(1 - p) * x + std::log(p);
+ }
+
+ // log of the binomial coefficient
+ static inline F log_binom_coeff(unsigned n, unsigned k) {
+ assert(n >= k);
+ if (n == k) return 0.0;
+ return lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1);
+ }
+
+ // http://en.wikipedia.org/wiki/Negative_binomial_distribution
+ // support [0, 1, 2 ...)
+ static inline F log_negative_binom(unsigned x, unsigned r, const F& p) {
+ assert(p > 0.0);
+ assert(p < 1.0);
+ return log_binom_coeff(x + r - 1u, x) + r * std::log(F(1) - p) + x * std::log(p);
+ }
+
+ // this is the Beta function, *not* the beta probability density
+ // http://mathworld.wolfram.com/BetaFunction.html
+ static inline F log_beta_fn(const F& x, const F& y) {
+ return lgamma(x) + lgamma(y) - lgamma(x + y);
+ }
+
+ // support x >= 0.0
+ static F log_gamma_density(const F& x, const F& shape, const F& rate) {
+ assert(x >= 0.0);
+ assert(shape > 0.0);
+ assert(rate > 0.0);
+ return (shape-1)*std::log(x) - shape*std::log(rate) - x/rate - lgamma(shape);
+ }
+
+ // this is the Beta *density* p(x ; alpha, beta)
+ // support x \in (0,1)
+ static inline F log_beta_density(const F& x, const F& alpha, const F& beta) {
+ assert(x > 0.0);
+ assert(x < 1.0);
+ assert(alpha > 0.0);
+ assert(beta > 0.0);
+ return (alpha-1)*std::log(x)+(beta-1)*std::log(1-x) - log_beta_fn(alpha, beta);
+ }
+
+ // note: this has been adapted so that 0 is in the support of the distribution
+ // support [0, 1, 2 ...)
+ static inline F log_yule_simon(unsigned x, const F& rho) {
+ assert(rho > 0.0);
+ return std::log(rho) + log_beta_fn(x + 1, rho + 1);
+ }
+
+ // see http://www.gatsby.ucl.ac.uk/~ywteh/research/compling/hpylm.pdf
+ // when y=1, sometimes written x^{\overline{n}} or x^{(n)} "Pochhammer symbol"
+ static inline F log_generalized_factorial(const F& x, const F& n, const F& y = 1.0) {
+ assert(x > 0.0);
+ assert(y >= 0.0);
+ assert(n > 0.0);
+ if (!n) return 0.0;
+ if (y == F(1)) {
+ return lgamma(x + n) - lgamma(x);
+ } else if (y) {
+ return n * std::log(y) + lgamma(x/y + n) - lgamma(x/y);
+ } else { // y == 0.0
+ return n * std::log(x);
+ }
+ }
+
+};
+
+typedef M<double> Md;
+typedef M<double> Mf;
+
+#endif
diff --git a/utils/m_test.cc b/utils/m_test.cc
new file mode 100644
index 00000000..fca8f895
--- /dev/null
+++ b/utils/m_test.cc
@@ -0,0 +1,75 @@
+#include "m.h"
+
+#include <iostream>
+#include <gtest/gtest.h>
+#include <cassert>
+
+using namespace std;
+
+class MTest : public testing::Test {
+ public:
+ MTest() {}
+ protected:
+ virtual void SetUp() { }
+ virtual void TearDown() { }
+};
+
+TEST_F(MTest, Poisson) {
+ double prev = 1.0;
+ double tot = 0;
+ for (int i = 0; i < 10; ++i) {
+ double p = Md::log_poisson(i, 0.99);
+ cerr << "p(i=" << i << ") = " << exp(p) << endl;
+ EXPECT_LT(p, prev);
+ tot += exp(p);
+ prev = p;
+ }
+ cerr << " tot=" << tot << endl;
+ EXPECT_LE(tot, 1.0);
+}
+
+TEST_F(MTest, YuleSimon) {
+ double prev = 1.0;
+ double tot = 0;
+ for (int i = 0; i < 10; ++i) {
+ double p = Md::log_yule_simon(i, 1.0);
+ cerr << "p(i=" << i << ") = " << exp(p) << endl;
+ EXPECT_LT(p, prev);
+ tot += exp(p);
+ prev = p;
+ }
+ cerr << " tot=" << tot << endl;
+ EXPECT_LE(tot, 1.0);
+}
+
+TEST_F(MTest, LogGeometric) {
+ double prev = 1.0;
+ double tot = 0;
+ for (int i = 0; i < 10; ++i) {
+ double p = Md::log_geometric(i, 0.5);
+ cerr << "p(i=" << i << ") = " << exp(p) << endl;
+ EXPECT_LT(p, prev);
+ tot += exp(p);
+ prev = p;
+ }
+ cerr << " tot=" << tot << endl;
+ EXPECT_LE(tot, 1.0);
+}
+
+TEST_F(MTest, GeneralizedFactorial) {
+ for (double i = 0.3; i < 10000; i += 0.4) {
+ double a = Md::log_generalized_factorial(1.0, i);
+ double b = lgamma(1.0 + i);
+ EXPECT_FLOAT_EQ(a,b);
+ }
+ double gf_3_6 = 3.0 * 4.0 * 5.0 * 6.0 * 7.0 * 8.0;
+ EXPECT_FLOAT_EQ(Md::log_generalized_factorial(3.0, 6.0), std::log(gf_3_6));
+ double gf_314_6 = 3.14 * 4.14 * 5.14 * 6.14 * 7.14 * 8.14;
+ EXPECT_FLOAT_EQ(Md::log_generalized_factorial(3.14, 6.0), std::log(gf_314_6));
+}
+
+int main(int argc, char** argv) {
+ testing::InitGoogleTest(&argc, argv);
+ return RUN_ALL_TESTS();
+}
+
diff --git a/utils/mfcr.h b/utils/mfcr.h
index 3eb133fc..396d0205 100644
--- a/utils/mfcr.h
+++ b/utils/mfcr.h
@@ -12,6 +12,7 @@
#include <boost/functional/hash.hpp>
#include "sampler.h"
#include "slice_sampler.h"
+#include "m.h"
struct TableCount {
TableCount() : count(), floor() {}
@@ -218,31 +219,14 @@ class MFCR {
return log_crp_prob(d_, alpha_);
}
- static double log_beta_density(const double& x, const double& alpha, const double& beta) {
- assert(x > 0.0);
- assert(x < 1.0);
- assert(alpha > 0.0);
- assert(beta > 0.0);
- const double lp = (alpha-1)*log(x)+(beta-1)*log(1-x)+lgamma(alpha+beta)-lgamma(alpha)-lgamma(beta);
- return lp;
- }
-
- static double log_gamma_density(const double& x, const double& shape, const double& rate) {
- assert(x >= 0.0);
- assert(shape > 0.0);
- assert(rate > 0.0);
- const double lp = (shape-1)*log(x) - shape*log(rate) - x/rate - lgamma(shape);
- return lp;
- }
-
// taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
// does not include draws from G_w's
double log_crp_prob(const double& d, const double& alpha) const {
double lp = 0.0;
if (has_d_prior())
- lp = log_beta_density(d, d_prior_alpha_, d_prior_beta_);
+ lp = Md::log_beta_density(d, d_prior_alpha_, d_prior_beta_);
if (has_alpha_prior())
- lp += log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
+ lp += Md::log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
assert(lp <= 0.0);
if (num_customers_) {
if (d > 0.0) {
--
cgit v1.2.3
From 8f6006cabee490a956940765c30cdd720d2e9161 Mon Sep 17 00:00:00 2001
From: Chris Dyer <cdyer@cs.cmu.edu>
Date: Sat, 3 Mar 2012 17:16:58 -0500
Subject: pyp lm, fixed hyperparameters inference
---
gi/pf/align-lexonly-pyp.cc | 2 +-
gi/pf/align-lexonly.cc | 2 +-
gi/pf/brat.cc | 2 +-
gi/pf/conditional_pseg.h | 4 +-
gi/pf/learn_cfg.cc | 4 +-
gi/pf/pfbrat.cc | 2 +-
gi/pf/pyp_lm.cc | 70 ++++++++++++++++++++++++++++---
phrasinator/gibbs_train_plm.cc | 2 +-
utils/ccrp.h | 95 ++++++++++++++++++------------------------
utils/ccrp_nt.h | 52 +++++++++++------------
utils/ccrp_onetable.h | 70 +++++++++++++++----------------
utils/mfcr.h | 58 +++++++++++++-------------
12 files changed, 203 insertions(+), 160 deletions(-)
(limited to 'utils/mfcr.h')
diff --git a/gi/pf/align-lexonly-pyp.cc b/gi/pf/align-lexonly-pyp.cc
index e24cb457..4ce7cf62 100644
--- a/gi/pf/align-lexonly-pyp.cc
+++ b/gi/pf/align-lexonly-pyp.cc
@@ -104,7 +104,7 @@ struct HierarchicalWordBase {
}
void Summary() const {
- cerr << "NUMBER OF CUSTOMERS: " << r.num_customers() << " (d=" << r.d() << ",\\alpha=" << r.alpha() << ')' << endl;
+ cerr << "NUMBER OF CUSTOMERS: " << r.num_customers() << " (d=" << r.discount() << ",\\alpha=" << r.alpha() << ')' << endl;
for (MFCR<vector<WordID> >::const_iterator it = r.begin(); it != r.end(); ++it)
cerr << " " << it->second.total_dish_count_ << " (on " << it->second.table_counts_.size() << " tables)" << TD::GetString(it->first) << endl;
}
diff --git a/gi/pf/align-lexonly.cc b/gi/pf/align-lexonly.cc
index 8c1d689f..dbc9dc07 100644
--- a/gi/pf/align-lexonly.cc
+++ b/gi/pf/align-lexonly.cc
@@ -105,7 +105,7 @@ struct HierarchicalWordBase {
}
void Summary() const {
- cerr << "NUMBER OF CUSTOMERS: " << r.num_customers() << " (\\alpha=" << r.concentration() << ')' << endl;
+ cerr << "NUMBER OF CUSTOMERS: " << r.num_customers() << " (\\alpha=" << r.alpha() << ')' << endl;
for (CCRP_NoTable<vector<WordID> >::const_iterator it = r.begin(); it != r.end(); ++it)
cerr << " " << it->second << '\t' << TD::GetString(it->first) << endl;
}
diff --git a/gi/pf/brat.cc b/gi/pf/brat.cc
index 7b60ef23..c2c52760 100644
--- a/gi/pf/brat.cc
+++ b/gi/pf/brat.cc
@@ -191,7 +191,7 @@ struct UniphraseLM {
void ResampleHyperparameters(MT19937* rng) {
phrases_.resample_hyperparameters(rng);
gen_.resample_hyperparameters(rng);
- cerr << " " << phrases_.concentration();
+ cerr << " " << phrases_.alpha();
}
CCRP_NoTable<vector<int> > phrases_;
diff --git a/gi/pf/conditional_pseg.h b/gi/pf/conditional_pseg.h
index 2e9e38fc..f9841cbf 100644
--- a/gi/pf/conditional_pseg.h
+++ b/gi/pf/conditional_pseg.h
@@ -22,7 +22,7 @@ struct MConditionalTranslationModel {
void Summary() const {
std::cerr << "Number of conditioning contexts: " << r.size() << std::endl;
for (RuleModelHash::const_iterator it = r.begin(); it != r.end(); ++it) {
- std::cerr << TD::GetString(it->first) << " \t(d=" << it->second.d() << ",\\alpha = " << it->second.alpha() << ") --------------------------" << std::endl;
+ std::cerr << TD::GetString(it->first) << " \t(d=" << it->second.discount() << ",\\alpha = " << it->second.alpha() << ") --------------------------" << std::endl;
for (MFCR<TRule>::const_iterator i2 = it->second.begin(); i2 != it->second.end(); ++i2)
std::cerr << " " << -1 << '\t' << i2->first << std::endl;
}
@@ -95,7 +95,7 @@ struct ConditionalTranslationModel {
void Summary() const {
std::cerr << "Number of conditioning contexts: " << r.size() << std::endl;
for (RuleModelHash::const_iterator it = r.begin(); it != r.end(); ++it) {
- std::cerr << TD::GetString(it->first) << " \t(\\alpha = " << it->second.concentration() << ") --------------------------" << std::endl;
+ std::cerr << TD::GetString(it->first) << " \t(\\alpha = " << it->second.alpha() << ") --------------------------" << std::endl;
for (CCRP_NoTable<TRule>::const_iterator i2 = it->second.begin(); i2 != it->second.end(); ++i2)
std::cerr << " " << i2->second << '\t' << i2->first << std::endl;
}
diff --git a/gi/pf/learn_cfg.cc b/gi/pf/learn_cfg.cc
index b2ca029a..5b748311 100644
--- a/gi/pf/learn_cfg.cc
+++ b/gi/pf/learn_cfg.cc
@@ -183,9 +183,9 @@ struct HieroLMModel {
nts[i].resample_hyperparameters(rng);
if (kHIERARCHICAL_PRIOR) {
q0.resample_hyperparameters(rng);
- cerr << "[base d=" << q0.discount() << ", alpha=" << q0.discount() << "]";
+ cerr << "[base d=" << q0.discount() << ", alpha=" << q0.alpha() << "]";
}
- cerr << " d=" << nts[0].discount() << ", alpha=" << nts[0].concentration() << endl;
+ cerr << " d=" << nts[0].discount() << ", alpha=" << nts[0].alpha() << endl;
}
const BaseRuleModel base;
diff --git a/gi/pf/pfbrat.cc b/gi/pf/pfbrat.cc
index 7b60ef23..c2c52760 100644
--- a/gi/pf/pfbrat.cc
+++ b/gi/pf/pfbrat.cc
@@ -191,7 +191,7 @@ struct UniphraseLM {
void ResampleHyperparameters(MT19937* rng) {
phrases_.resample_hyperparameters(rng);
gen_.resample_hyperparameters(rng);
- cerr << " " << phrases_.concentration();
+ cerr << " " << phrases_.alpha();
}
CCRP_NoTable<vector<int> > phrases_;
diff --git a/gi/pf/pyp_lm.cc b/gi/pf/pyp_lm.cc
index 2837e33c..0d85536c 100644
--- a/gi/pf/pyp_lm.cc
+++ b/gi/pf/pyp_lm.cc
@@ -50,16 +50,19 @@ template <unsigned N> struct PYPLM;
// uniform base distribution
template<> struct PYPLM<0> {
- PYPLM(unsigned vs) : p0(1.0 / vs) {}
- void increment(WordID w, const vector<WordID>& context, MT19937* rng) const {}
- void decrement(WordID w, const vector<WordID>& context, MT19937* rng) const {}
+ PYPLM(unsigned vs) : p0(1.0 / vs), draws() {}
+ void increment(WordID w, const vector<WordID>& context, MT19937* rng) { ++draws; }
+ void decrement(WordID w, const vector<WordID>& context, MT19937* rng) { --draws; assert(draws >= 0); }
double prob(WordID w, const vector<WordID>& context) const { return p0; }
+ void resample_hyperparameters(MT19937* rng, const unsigned nloop, const unsigned niterations) {}
+ double log_likelihood() const { return draws * log(p0); }
const double p0;
+ int draws;
};
// represents an N-gram LM
template <unsigned N> struct PYPLM {
- PYPLM(unsigned vs) : backoff(vs) {}
+ PYPLM(unsigned vs) : backoff(vs), d(0.8), alpha(1.0) {}
void increment(WordID w, const vector<WordID>& context, MT19937* rng) {
const double bo = backoff.prob(w, context);
static vector<WordID> lookup(N-1);
@@ -67,7 +70,7 @@ template <unsigned N> struct PYPLM {
lookup[i] = context[context.size() - 1 - i];
typename unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > >::iterator it = p.find(lookup);
if (it == p.end())
- it = p.insert(make_pair(lookup, CCRP<WordID>(1,1,1,1))).first;
+ it = p.insert(make_pair(lookup, CCRP<WordID>(d,alpha))).first;
if (it->second.increment(w, bo, rng))
backoff.increment(w, context, rng);
}
@@ -89,7 +92,58 @@ template <unsigned N> struct PYPLM {
if (it == p.end()) return bo;
return it->second.prob(w, bo);
}
+
+ double log_likelihood() const {
+ return log_likelihood(d, alpha) + backoff.log_likelihood();
+ }
+
+ double log_likelihood(const double& dd, const double& aa) const {
+ if (aa <= -dd) return -std::numeric_limits<double>::infinity();
+ double llh = Md::log_beta_density(dd, 1, 1) + Md::log_gamma_density(aa, 1, 1);
+ typename unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > >::const_iterator it;
+ for (it = p.begin(); it != p.end(); ++it)
+ llh += it->second.log_crp_prob(dd, aa);
+ return llh;
+ }
+
+ struct DiscountResampler {
+ DiscountResampler(const PYPLM& m) : m_(m) {}
+ const PYPLM& m_;
+ double operator()(const double& proposed_discount) const {
+ return m_.log_likelihood(proposed_discount, m_.alpha);
+ }
+ };
+
+ struct AlphaResampler {
+ AlphaResampler(const PYPLM& m) : m_(m) {}
+ const PYPLM& m_;
+ double operator()(const double& proposed_alpha) const {
+ return m_.log_likelihood(m_.d, proposed_alpha);
+ }
+ };
+
+ void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
+ DiscountResampler dr(*this);
+ AlphaResampler ar(*this);
+ for (int iter = 0; iter < nloop; ++iter) {
+ alpha = slice_sampler1d(ar, alpha, *rng, 0.0,
+ std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+ d = slice_sampler1d(dr, d, *rng, std::numeric_limits<double>::min(),
+ 1.0, 0.0, niterations, 100*niterations);
+ }
+ alpha = slice_sampler1d(ar, alpha, *rng, 0.0,
+ std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+ typename unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > >::iterator it;
+ cerr << "PYPLM<" << N << ">(d=" << d << ",a=" << alpha << ") = " << log_likelihood(d, alpha) << endl;
+ for (it = p.begin(); it != p.end(); ++it) {
+ it->second.set_discount(d);
+ it->second.set_alpha(alpha);
+ }
+ backoff.resample_hyperparameters(rng, nloop, niterations);
+ }
+
PYPLM<N-1> backoff;
+ double d, alpha;
unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > > p;
};
@@ -109,7 +163,7 @@ int main(int argc, char** argv) {
cerr << "Reading corpus...\n";
CorpusTools::ReadFromFile(conf["input"].as<string>(), &corpuse, &vocabe);
cerr << "E-corpus size: " << corpuse.size() << " sentences\t (" << vocabe.size() << " word types)\n";
-#define kORDER 5
+#define kORDER 3
PYPLM<kORDER> lm(vocabe.size());
vector<WordID> ctx(kORDER - 1, TD::Convert("<s>"));
int mci = corpuse.size() * 99 / 100;
@@ -126,6 +180,10 @@ int main(int argc, char** argv) {
if (SS > 0) lm.decrement(kEOS, ctx, &rng);
lm.increment(kEOS, ctx, &rng);
}
+ if (SS % 10 == 9) {
+ cerr << " [LLH=" << lm.log_likelihood() << "]" << endl;
+ if (SS % 20 == 19) lm.resample_hyperparameters(&rng);
+ } else { cerr << '.' << flush; }
}
double llh = 0;
unsigned cnt = 0;
diff --git a/phrasinator/gibbs_train_plm.cc b/phrasinator/gibbs_train_plm.cc
index 66b46011..54861dcb 100644
--- a/phrasinator/gibbs_train_plm.cc
+++ b/phrasinator/gibbs_train_plm.cc
@@ -252,7 +252,7 @@ struct UniphraseLM {
void ResampleHyperparameters(MT19937* rng) {
phrases_.resample_hyperparameters(rng);
gen_.resample_hyperparameters(rng);
- cerr << " d=" << phrases_.discount() << ",c=" << phrases_.concentration();
+ cerr << " d=" << phrases_.discount() << ",a=" << phrases_.alpha();
}
CCRP<vector<int> > phrases_;
diff --git a/utils/ccrp.h b/utils/ccrp.h
index 1a9e3ed5..d9a38089 100644
--- a/utils/ccrp.h
+++ b/utils/ccrp.h
@@ -17,35 +17,37 @@
template <typename Dish, typename DishHash = boost::hash<Dish> >
class CCRP {
public:
- CCRP(double disc, double conc) :
+ CCRP(double disc, double alpha) :
num_tables_(),
num_customers_(),
discount_(disc),
- concentration_(conc),
+ alpha_(alpha),
discount_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
- concentration_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
- concentration_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+ alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+ alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
CCRP(double d_alpha, double d_beta, double c_shape, double c_rate, double d = 0.9, double c = 1.0) :
num_tables_(),
num_customers_(),
discount_(d),
- concentration_(c),
+ alpha_(c),
discount_prior_alpha_(d_alpha),
discount_prior_beta_(d_beta),
- concentration_prior_shape_(c_shape),
- concentration_prior_rate_(c_rate) {}
+ alpha_prior_shape_(c_shape),
+ alpha_prior_rate_(c_rate) {}
double discount() const { return discount_; }
- double concentration() const { return concentration_; }
+ double alpha() const { return alpha_; }
+ void set_discount(double d) { discount_ = d; }
+ void set_alpha(double a) { alpha_ = a; }
bool has_discount_prior() const {
return !std::isnan(discount_prior_alpha_);
}
- bool has_concentration_prior() const {
- return !std::isnan(concentration_prior_shape_);
+ bool has_alpha_prior() const {
+ return !std::isnan(alpha_prior_shape_);
}
void clear() {
@@ -79,7 +81,7 @@ class CCRP {
DishLocations& loc = dish_locs_[dish];
bool share_table = false;
if (loc.total_dish_count_) {
- const double p_empty = (concentration_ + num_tables_ * discount_) * p0;
+ const double p_empty = (alpha_ + num_tables_ * discount_) * p0;
const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
share_table = rng->SelectSample(p_empty, p_share);
}
@@ -113,7 +115,7 @@ class CCRP {
DishLocations& loc = dish_locs_[dish];
bool share_table = false;
if (loc.total_dish_count_) {
- const T p_empty = T(concentration_ + num_tables_ * discount_) * p0;
+ const T p_empty = T(alpha_ + num_tables_ * discount_) * p0;
const T p_share = T(loc.total_dish_count_ - loc.table_counts_.size() * discount_);
share_table = rng->SelectSample(p_empty, p_share);
}
@@ -180,63 +182,46 @@ class CCRP {
double prob(const Dish& dish, const double& p0) const {
const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
- const double r = num_tables_ * discount_ + concentration_;
+ const double r = num_tables_ * discount_ + alpha_;
if (it == dish_locs_.end()) {
- return r * p0 / (num_customers_ + concentration_);
+ return r * p0 / (num_customers_ + alpha_);
} else {
return (it->second.total_dish_count_ - discount_ * it->second.table_counts_.size() + r * p0) /
- (num_customers_ + concentration_);
+ (num_customers_ + alpha_);
}
}
template <typename T>
T probT(const Dish& dish, const T& p0) const {
const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
- const T r = T(num_tables_ * discount_ + concentration_);
+ const T r = T(num_tables_ * discount_ + alpha_);
if (it == dish_locs_.end()) {
- return r * p0 / T(num_customers_ + concentration_);
+ return r * p0 / T(num_customers_ + alpha_);
} else {
return (T(it->second.total_dish_count_ - discount_ * it->second.table_counts_.size()) + r * p0) /
- T(num_customers_ + concentration_);
+ T(num_customers_ + alpha_);
}
}
double log_crp_prob() const {
- return log_crp_prob(discount_, concentration_);
- }
-
- static double log_beta_density(const double& x, const double& alpha, const double& beta) {
- assert(x > 0.0);
- assert(x < 1.0);
- assert(alpha > 0.0);
- assert(beta > 0.0);
- const double lp = (alpha-1)*log(x)+(beta-1)*log(1-x)+lgamma(alpha+beta)-lgamma(alpha)-lgamma(beta);
- return lp;
- }
-
- static double log_gamma_density(const double& x, const double& shape, const double& rate) {
- assert(x >= 0.0);
- assert(shape > 0.0);
- assert(rate > 0.0);
- const double lp = (shape-1)*log(x) - shape*log(rate) - x/rate - lgamma(shape);
- return lp;
+ return log_crp_prob(discount_, alpha_);
}
// taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
// does not include P_0's
- double log_crp_prob(const double& discount, const double& concentration) const {
+ double log_crp_prob(const double& discount, const double& alpha) const {
double lp = 0.0;
if (has_discount_prior())
- lp = log_beta_density(discount, discount_prior_alpha_, discount_prior_beta_);
- if (has_concentration_prior())
- lp += log_gamma_density(concentration, concentration_prior_shape_, concentration_prior_rate_);
+ lp = Md::log_beta_density(discount, discount_prior_alpha_, discount_prior_beta_);
+ if (has_alpha_prior())
+ lp += Md::log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
assert(lp <= 0.0);
if (num_customers_) {
if (discount > 0.0) {
const double r = lgamma(1.0 - discount);
- lp += lgamma(concentration) - lgamma(concentration + num_customers_)
- + num_tables_ * log(discount) + lgamma(concentration / discount + num_tables_)
- - lgamma(concentration / discount);
+ lp += lgamma(alpha) - lgamma(alpha + num_customers_)
+ + num_tables_ * log(discount) + lgamma(alpha / discount + num_tables_)
+ - lgamma(alpha / discount);
assert(std::isfinite(lp));
for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
it != dish_locs_.end(); ++it) {
@@ -254,12 +239,12 @@ class CCRP {
}
void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
- assert(has_discount_prior() || has_concentration_prior());
+ assert(has_discount_prior() || has_alpha_prior());
DiscountResampler dr(*this);
ConcentrationResampler cr(*this);
for (int iter = 0; iter < nloop; ++iter) {
- if (has_concentration_prior()) {
- concentration_ = slice_sampler1d(cr, concentration_, *rng, 0.0,
+ if (has_alpha_prior()) {
+ alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
if (has_discount_prior()) {
@@ -267,7 +252,7 @@ class CCRP {
1.0, 0.0, niterations, 100*niterations);
}
}
- concentration_ = slice_sampler1d(cr, concentration_, *rng, 0.0,
+ alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
@@ -275,15 +260,15 @@ class CCRP {
DiscountResampler(const CCRP& crp) : crp_(crp) {}
const CCRP& crp_;
double operator()(const double& proposed_discount) const {
- return crp_.log_crp_prob(proposed_discount, crp_.concentration_);
+ return crp_.log_crp_prob(proposed_discount, crp_.alpha_);
}
};
struct ConcentrationResampler {
ConcentrationResampler(const CCRP& crp) : crp_(crp) {}
const CCRP& crp_;
- double operator()(const double& proposed_concentration) const {
- return crp_.log_crp_prob(crp_.discount_, proposed_concentration);
+ double operator()(const double& proposed_alpha) const {
+ return crp_.log_crp_prob(crp_.discount_, proposed_alpha);
}
};
@@ -295,7 +280,7 @@ class CCRP {
};
void Print(std::ostream* out) const {
- std::cerr << "PYP(d=" << discount_ << ",c=" << concentration_ << ") customers=" << num_customers_ << std::endl;
+ std::cerr << "PYP(d=" << discount_ << ",c=" << alpha_ << ") customers=" << num_customers_ << std::endl;
for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
it != dish_locs_.end(); ++it) {
(*out) << it->first << " (" << it->second.total_dish_count_ << " on " << it->second.table_counts_.size() << " tables): ";
@@ -320,15 +305,15 @@ class CCRP {
std::tr1::unordered_map<Dish, DishLocations, DishHash> dish_locs_;
double discount_;
- double concentration_;
+ double alpha_;
// optional beta prior on discount_ (NaN if no prior)
double discount_prior_alpha_;
double discount_prior_beta_;
- // optional gamma prior on concentration_ (NaN if no prior)
- double concentration_prior_shape_;
- double concentration_prior_rate_;
+ // optional gamma prior on alpha_ (NaN if no prior)
+ double alpha_prior_shape_;
+ double alpha_prior_rate_;
};
template <typename T,typename H>
diff --git a/utils/ccrp_nt.h b/utils/ccrp_nt.h
index 63b6f4c2..79321493 100644
--- a/utils/ccrp_nt.h
+++ b/utils/ccrp_nt.h
@@ -18,20 +18,20 @@ class CCRP_NoTable {
public:
explicit CCRP_NoTable(double conc) :
num_customers_(),
- concentration_(conc),
- concentration_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
- concentration_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+ alpha_(conc),
+ alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+ alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
CCRP_NoTable(double c_shape, double c_rate, double c = 10.0) :
num_customers_(),
- concentration_(c),
- concentration_prior_shape_(c_shape),
- concentration_prior_rate_(c_rate) {}
+ alpha_(c),
+ alpha_prior_shape_(c_shape),
+ alpha_prior_rate_(c_rate) {}
- double concentration() const { return concentration_; }
+ double alpha() const { return alpha_; }
- bool has_concentration_prior() const {
- return !std::isnan(concentration_prior_shape_);
+ bool has_alpha_prior() const {
+ return !std::isnan(alpha_prior_shape_);
}
void clear() {
@@ -73,16 +73,16 @@ class CCRP_NoTable {
double prob(const Dish& dish, const double& p0) const {
const unsigned at_table = num_customers(dish);
- return (at_table + p0 * concentration_) / (num_customers_ + concentration_);
+ return (at_table + p0 * alpha_) / (num_customers_ + alpha_);
}
double logprob(const Dish& dish, const double& logp0) const {
const unsigned at_table = num_customers(dish);
- return log(at_table + exp(logp0 + log(concentration_))) - log(num_customers_ + concentration_);
+ return log(at_table + exp(logp0 + log(alpha_))) - log(num_customers_ + alpha_);
}
double log_crp_prob() const {
- return log_crp_prob(concentration_);
+ return log_crp_prob(alpha_);
}
static double log_gamma_density(const double& x, const double& shape, const double& rate) {
@@ -95,14 +95,14 @@ class CCRP_NoTable {
// taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
// does not include P_0's
- double log_crp_prob(const double& concentration) const {
+ double log_crp_prob(const double& alpha) const {
double lp = 0.0;
- if (has_concentration_prior())
- lp += log_gamma_density(concentration, concentration_prior_shape_, concentration_prior_rate_);
+ if (has_alpha_prior())
+ lp += log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
assert(lp <= 0.0);
if (num_customers_) {
- lp += lgamma(concentration) - lgamma(concentration + num_customers_) +
- custs_.size() * log(concentration);
+ lp += lgamma(alpha) - lgamma(alpha + num_customers_) +
+ custs_.size() * log(alpha);
assert(std::isfinite(lp));
for (typename std::tr1::unordered_map<Dish, unsigned, DishHash>::const_iterator it = custs_.begin();
it != custs_.end(); ++it) {
@@ -114,10 +114,10 @@ class CCRP_NoTable {
}
void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
- assert(has_concentration_prior());
+ assert(has_alpha_prior());
ConcentrationResampler cr(*this);
for (int iter = 0; iter < nloop; ++iter) {
- concentration_ = slice_sampler1d(cr, concentration_, *rng, 0.0,
+ alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
}
@@ -125,13 +125,13 @@ class CCRP_NoTable {
struct ConcentrationResampler {
ConcentrationResampler(const CCRP_NoTable& crp) : crp_(crp) {}
const CCRP_NoTable& crp_;
- double operator()(const double& proposed_concentration) const {
- return crp_.log_crp_prob(proposed_concentration);
+ double operator()(const double& proposed_alpha) const {
+ return crp_.log_crp_prob(proposed_alpha);
}
};
void Print(std::ostream* out) const {
- (*out) << "DP(alpha=" << concentration_ << ") customers=" << num_customers_ << std::endl;
+ (*out) << "DP(alpha=" << alpha_ << ") customers=" << num_customers_ << std::endl;
int cc = 0;
for (typename std::tr1::unordered_map<Dish, unsigned, DishHash>::const_iterator it = custs_.begin();
it != custs_.end(); ++it) {
@@ -153,11 +153,11 @@ class CCRP_NoTable {
return custs_.end();
}
- double concentration_;
+ double alpha_;
- // optional gamma prior on concentration_ (NaN if no prior)
- double concentration_prior_shape_;
- double concentration_prior_rate_;
+ // optional gamma prior on alpha_ (NaN if no prior)
+ double alpha_prior_shape_;
+ double alpha_prior_rate_;
};
template <typename T,typename H>
diff --git a/utils/ccrp_onetable.h b/utils/ccrp_onetable.h
index b63737d1..1fe01b0e 100644
--- a/utils/ccrp_onetable.h
+++ b/utils/ccrp_onetable.h
@@ -21,33 +21,33 @@ class CCRP_OneTable {
num_tables_(),
num_customers_(),
discount_(disc),
- concentration_(conc),
+ alpha_(conc),
discount_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
- concentration_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
- concentration_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+ alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+ alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
CCRP_OneTable(double d_alpha, double d_beta, double c_shape, double c_rate, double d = 0.9, double c = 1.0) :
num_tables_(),
num_customers_(),
discount_(d),
- concentration_(c),
+ alpha_(c),
discount_prior_alpha_(d_alpha),
discount_prior_beta_(d_beta),
- concentration_prior_shape_(c_shape),
- concentration_prior_rate_(c_rate) {}
+ alpha_prior_shape_(c_shape),
+ alpha_prior_rate_(c_rate) {}
double discount() const { return discount_; }
- double concentration() const { return concentration_; }
- void set_concentration(double c) { concentration_ = c; }
+ double alpha() const { return alpha_; }
+ void set_alpha(double c) { alpha_ = c; }
void set_discount(double d) { discount_ = d; }
bool has_discount_prior() const {
return !std::isnan(discount_prior_alpha_);
}
- bool has_concentration_prior() const {
- return !std::isnan(concentration_prior_shape_);
+ bool has_alpha_prior() const {
+ return !std::isnan(alpha_prior_shape_);
}
void clear() {
@@ -108,29 +108,29 @@ class CCRP_OneTable {
double prob(const Dish& dish, const double& p0) const {
const typename DishMapType::const_iterator it = dish_counts_.find(dish);
- const double r = num_tables_ * discount_ + concentration_;
+ const double r = num_tables_ * discount_ + alpha_;
if (it == dish_counts_.end()) {
- return r * p0 / (num_customers_ + concentration_);
+ return r * p0 / (num_customers_ + alpha_);
} else {
return (it->second - discount_ + r * p0) /
- (num_customers_ + concentration_);
+ (num_customers_ + alpha_);
}
}
template <typename T>
T probT(const Dish& dish, const T& p0) const {
const typename DishMapType::const_iterator it = dish_counts_.find(dish);
- const T r(num_tables_ * discount_ + concentration_);
+ const T r(num_tables_ * discount_ + alpha_);
if (it == dish_counts_.end()) {
- return r * p0 / T(num_customers_ + concentration_);
+ return r * p0 / T(num_customers_ + alpha_);
} else {
return (T(it->second - discount_) + r * p0) /
- T(num_customers_ + concentration_);
+ T(num_customers_ + alpha_);
}
}
double log_crp_prob() const {
- return log_crp_prob(discount_, concentration_);
+ return log_crp_prob(discount_, alpha_);
}
static double log_beta_density(const double& x, const double& alpha, const double& beta) {
@@ -152,19 +152,19 @@ class CCRP_OneTable {
// taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
// does not include P_0's
- double log_crp_prob(const double& discount, const double& concentration) const {
+ double log_crp_prob(const double& discount, const double& alpha) const {
double lp = 0.0;
if (has_discount_prior())
lp = log_beta_density(discount, discount_prior_alpha_, discount_prior_beta_);
- if (has_concentration_prior())
- lp += log_gamma_density(concentration, concentration_prior_shape_, concentration_prior_rate_);
+ if (has_alpha_prior())
+ lp += log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
assert(lp <= 0.0);
if (num_customers_) {
if (discount > 0.0) {
const double r = lgamma(1.0 - discount);
- lp += lgamma(concentration) - lgamma(concentration + num_customers_)
- + num_tables_ * log(discount) + lgamma(concentration / discount + num_tables_)
- - lgamma(concentration / discount);
+ lp += lgamma(alpha) - lgamma(alpha + num_customers_)
+ + num_tables_ * log(discount) + lgamma(alpha / discount + num_tables_)
+ - lgamma(alpha / discount);
assert(std::isfinite(lp));
for (typename DishMapType::const_iterator it = dish_counts_.begin();
it != dish_counts_.end(); ++it) {
@@ -180,12 +180,12 @@ class CCRP_OneTable {
}
void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
- assert(has_discount_prior() || has_concentration_prior());
+ assert(has_discount_prior() || has_alpha_prior());
DiscountResampler dr(*this);
ConcentrationResampler cr(*this);
for (int iter = 0; iter < nloop; ++iter) {
- if (has_concentration_prior()) {
- concentration_ = slice_sampler1d(cr, concentration_, *rng, 0.0,
+ if (has_alpha_prior()) {
+ alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
if (has_discount_prior()) {
@@ -193,7 +193,7 @@ class CCRP_OneTable {
1.0, 0.0, niterations, 100*niterations);
}
}
- concentration_ = slice_sampler1d(cr, concentration_, *rng, 0.0,
+ alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
@@ -201,20 +201,20 @@ class CCRP_OneTable {
DiscountResampler(const CCRP_OneTable& crp) : crp_(crp) {}
const CCRP_OneTable& crp_;
double operator()(const double& proposed_discount) const {
- return crp_.log_crp_prob(proposed_discount, crp_.concentration_);
+ return crp_.log_crp_prob(proposed_discount, crp_.alpha_);
}
};
struct ConcentrationResampler {
ConcentrationResampler(const CCRP_OneTable& crp) : crp_(crp) {}
const CCRP_OneTable& crp_;
- double operator()(const double& proposed_concentration) const {
- return crp_.log_crp_prob(crp_.discount_, proposed_concentration);
+ double operator()(const double& proposed_alpha) const {
+ return crp_.log_crp_prob(crp_.discount_, proposed_alpha);
}
};
void Print(std::ostream* out) const {
- (*out) << "PYP(d=" << discount_ << ",c=" << concentration_ << ") customers=" << num_customers_ << std::endl;
+ (*out) << "PYP(d=" << discount_ << ",c=" << alpha_ << ") customers=" << num_customers_ << std::endl;
for (typename DishMapType::const_iterator it = dish_counts_.begin(); it != dish_counts_.end(); ++it) {
(*out) << " " << it->first << " = " << it->second << std::endl;
}
@@ -233,15 +233,15 @@ class CCRP_OneTable {
DishMapType dish_counts_;
double discount_;
- double concentration_;
+ double alpha_;
// optional beta prior on discount_ (NaN if no prior)
double discount_prior_alpha_;
double discount_prior_beta_;
- // optional gamma prior on concentration_ (NaN if no prior)
- double concentration_prior_shape_;
- double concentration_prior_rate_;
+ // optional gamma prior on alpha_ (NaN if no prior)
+ double alpha_prior_shape_;
+ double alpha_prior_rate_;
};
template <typename T,typename H>
diff --git a/utils/mfcr.h b/utils/mfcr.h
index 396d0205..df988f51 100644
--- a/utils/mfcr.h
+++ b/utils/mfcr.h
@@ -43,29 +43,29 @@ class MFCR {
num_floors_(num_floors),
num_tables_(),
num_customers_(),
- d_(d),
+ discount_(d),
alpha_(alpha),
- d_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
- d_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
+ discount_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
+ discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
- MFCR(unsigned num_floors, double d_alpha, double d_beta, double alpha_shape, double alpha_rate, double d = 0.9, double alpha = 10.0) :
+ MFCR(unsigned num_floors, double discount_alpha, double discount_beta, double alpha_shape, double alpha_rate, double d = 0.9, double alpha = 10.0) :
num_floors_(num_floors),
num_tables_(),
num_customers_(),
- d_(d),
+ discount_(d),
alpha_(alpha),
- d_prior_alpha_(d_alpha),
- d_prior_beta_(d_beta),
+ discount_prior_alpha_(discount_alpha),
+ discount_prior_beta_(discount_beta),
alpha_prior_shape_(alpha_shape),
alpha_prior_rate_(alpha_rate) {}
- double d() const { return d_; }
+ double discount() const { return discount_; }
double alpha() const { return alpha_; }
- bool has_d_prior() const {
- return !std::isnan(d_prior_alpha_);
+ bool has_discount_prior() const {
+ return !std::isnan(discount_prior_alpha_);
}
bool has_alpha_prior() const {
@@ -122,15 +122,15 @@ class MFCR {
int floor = -1;
bool share_table = false;
if (loc.total_dish_count_) {
- const double p_empty = (alpha_ + num_tables_ * d_) * marg_p0;
- const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * d_);
+ const double p_empty = (alpha_ + num_tables_ * discount_) * marg_p0;
+ const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
share_table = rng->SelectSample(p_empty, p_share);
}
if (share_table) {
- double r = rng->next() * (loc.total_dish_count_ - loc.table_counts_.size() * d_);
+ double r = rng->next() * (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
for (typename std::list<TableCount>::iterator ti = loc.table_counts_.begin();
ti != loc.table_counts_.end(); ++ti) {
- r -= ti->count - d_;
+ r -= ti->count - discount_;
if (r <= 0.0) {
++ti->count;
floor = ti->floor;
@@ -206,25 +206,25 @@ class MFCR {
const double marg_p0 = std::inner_product(p0s.begin(), p0s.end(), lambdas.begin(), 0.0);
assert(marg_p0 <= 1.0);
const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
- const double r = num_tables_ * d_ + alpha_;
+ const double r = num_tables_ * discount_ + alpha_;
if (it == dish_locs_.end()) {
return r * marg_p0 / (num_customers_ + alpha_);
} else {
- return (it->second.total_dish_count_ - d_ * it->second.table_counts_.size() + r * marg_p0) /
+ return (it->second.total_dish_count_ - discount_ * it->second.table_counts_.size() + r * marg_p0) /
(num_customers_ + alpha_);
}
}
double log_crp_prob() const {
- return log_crp_prob(d_, alpha_);
+ return log_crp_prob(discount_, alpha_);
}
// taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
// does not include draws from G_w's
double log_crp_prob(const double& d, const double& alpha) const {
double lp = 0.0;
- if (has_d_prior())
- lp = Md::log_beta_density(d, d_prior_alpha_, d_prior_beta_);
+ if (has_discount_prior())
+ lp = Md::log_beta_density(d, discount_prior_alpha_, discount_prior_beta_);
if (has_alpha_prior())
lp += Md::log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
assert(lp <= 0.0);
@@ -251,7 +251,7 @@ class MFCR {
}
void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
- assert(has_d_prior() || has_alpha_prior());
+ assert(has_discount_prior() || has_alpha_prior());
DiscountResampler dr(*this);
ConcentrationResampler cr(*this);
for (int iter = 0; iter < nloop; ++iter) {
@@ -259,8 +259,8 @@ class MFCR {
alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
- if (has_d_prior()) {
- d_ = slice_sampler1d(dr, d_, *rng, std::numeric_limits<double>::min(),
+ if (has_discount_prior()) {
+ discount_ = slice_sampler1d(dr, discount_, *rng, std::numeric_limits<double>::min(),
1.0, 0.0, niterations, 100*niterations);
}
}
@@ -279,8 +279,8 @@ class MFCR {
struct ConcentrationResampler {
ConcentrationResampler(const MFCR& crp) : crp_(crp) {}
const MFCR& crp_;
- double operator()(const double& proposed_alpha) const {
- return crp_.log_crp_prob(crp_.d_, proposed_alpha);
+ double operator()(const double& proposediscount_alpha) const {
+ return crp_.log_crp_prob(crp_.discount_, proposediscount_alpha);
}
};
@@ -292,7 +292,7 @@ class MFCR {
};
void Print(std::ostream* out) const {
- (*out) << "MFCR(d=" << d_ << ",alpha=" << alpha_ << ") customers=" << num_customers_ << std::endl;
+ (*out) << "MFCR(d=" << discount_ << ",alpha=" << alpha_ << ") customers=" << num_customers_ << std::endl;
for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
it != dish_locs_.end(); ++it) {
(*out) << it->first << " (" << it->second.total_dish_count_ << " on " << it->second.table_counts_.size() << " tables): ";
@@ -317,12 +317,12 @@ class MFCR {
unsigned num_customers_;
std::tr1::unordered_map<Dish, DishLocations, DishHash> dish_locs_;
- double d_;
+ double discount_;
double alpha_;
- // optional beta prior on d_ (NaN if no prior)
- double d_prior_alpha_;
- double d_prior_beta_;
+ // optional beta prior on discount_ (NaN if no prior)
+ double discount_prior_alpha_;
+ double discount_prior_beta_;
// optional gamma prior on alpha_ (NaN if no prior)
double alpha_prior_shape_;
--
cgit v1.2.3
From 1d5a0055a948663d799b4c5b1380ce1d9742bf6b Mon Sep 17 00:00:00 2001
From: Chris Dyer <cdyer@cs.cmu.edu>
Date: Mon, 5 Mar 2012 14:51:04 -0500
Subject: support strength=0 PYPs, final notation clean-up
---
gi/pf/align-lexonly-pyp.cc | 2 +-
gi/pf/conditional_pseg.h | 2 +-
gi/pf/learn_cfg.cc | 4 +-
gi/pf/pyp_lm.cc | 22 ++++-----
phrasinator/gibbs_train_plm.cc | 2 +-
utils/ccrp.h | 106 ++++++++++++++++++++++-------------------
utils/mfcr.h | 105 ++++++++++++++++++++++------------------
7 files changed, 131 insertions(+), 112 deletions(-)
(limited to 'utils/mfcr.h')
diff --git a/gi/pf/align-lexonly-pyp.cc b/gi/pf/align-lexonly-pyp.cc
index 4ce7cf62..87f7f6b5 100644
--- a/gi/pf/align-lexonly-pyp.cc
+++ b/gi/pf/align-lexonly-pyp.cc
@@ -104,7 +104,7 @@ struct HierarchicalWordBase {
}
void Summary() const {
- cerr << "NUMBER OF CUSTOMERS: " << r.num_customers() << " (d=" << r.discount() << ",\\alpha=" << r.alpha() << ')' << endl;
+ cerr << "NUMBER OF CUSTOMERS: " << r.num_customers() << " (d=" << r.discount() << ",s=" << r.strength() << ')' << endl;
for (MFCR<vector<WordID> >::const_iterator it = r.begin(); it != r.end(); ++it)
cerr << " " << it->second.total_dish_count_ << " (on " << it->second.table_counts_.size() << " tables)" << TD::GetString(it->first) << endl;
}
diff --git a/gi/pf/conditional_pseg.h b/gi/pf/conditional_pseg.h
index f9841cbf..86403d8d 100644
--- a/gi/pf/conditional_pseg.h
+++ b/gi/pf/conditional_pseg.h
@@ -22,7 +22,7 @@ struct MConditionalTranslationModel {
void Summary() const {
std::cerr << "Number of conditioning contexts: " << r.size() << std::endl;
for (RuleModelHash::const_iterator it = r.begin(); it != r.end(); ++it) {
- std::cerr << TD::GetString(it->first) << " \t(d=" << it->second.discount() << ",\\alpha = " << it->second.alpha() << ") --------------------------" << std::endl;
+ std::cerr << TD::GetString(it->first) << " \t(d=" << it->second.discount() << ",s=" << it->second.strength() << ") --------------------------" << std::endl;
for (MFCR<TRule>::const_iterator i2 = it->second.begin(); i2 != it->second.end(); ++i2)
std::cerr << " " << -1 << '\t' << i2->first << std::endl;
}
diff --git a/gi/pf/learn_cfg.cc b/gi/pf/learn_cfg.cc
index 5b748311..bf157828 100644
--- a/gi/pf/learn_cfg.cc
+++ b/gi/pf/learn_cfg.cc
@@ -183,9 +183,9 @@ struct HieroLMModel {
nts[i].resample_hyperparameters(rng);
if (kHIERARCHICAL_PRIOR) {
q0.resample_hyperparameters(rng);
- cerr << "[base d=" << q0.discount() << ", alpha=" << q0.alpha() << "]";
+ cerr << "[base d=" << q0.discount() << ", s=" << q0.strength() << "]";
}
- cerr << " d=" << nts[0].discount() << ", alpha=" << nts[0].alpha() << endl;
+ cerr << " d=" << nts[0].discount() << ", s=" << nts[0].strength() << endl;
}
const BaseRuleModel base;
diff --git a/gi/pf/pyp_lm.cc b/gi/pf/pyp_lm.cc
index e5c44c8b..7ebada13 100644
--- a/gi/pf/pyp_lm.cc
+++ b/gi/pf/pyp_lm.cc
@@ -78,14 +78,14 @@ template <unsigned N> struct PYPLM {
backoff(vs, da, db, ss, sr),
discount_a(da), discount_b(db),
strength_s(ss), strength_r(sr),
- d(0.8), alpha(1.0), lookup(N-1) {}
+ d(0.8), strength(1.0), lookup(N-1) {}
void increment(WordID w, const vector<WordID>& context, MT19937* rng) {
const double bo = backoff.prob(w, context);
for (unsigned i = 0; i < N-1; ++i)
lookup[i] = context[context.size() - 1 - i];
typename unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > >::iterator it = p.find(lookup);
if (it == p.end())
- it = p.insert(make_pair(lookup, CCRP<WordID>(d,alpha))).first;
+ it = p.insert(make_pair(lookup, CCRP<WordID>(d,strength))).first;
if (it->second.increment(w, bo, rng))
backoff.increment(w, context, rng);
}
@@ -107,7 +107,7 @@ template <unsigned N> struct PYPLM {
}
double log_likelihood() const {
- return log_likelihood(d, alpha) + backoff.log_likelihood();
+ return log_likelihood(d, strength) + backoff.log_likelihood();
}
double log_likelihood(const double& dd, const double& aa) const {
@@ -125,15 +125,15 @@ template <unsigned N> struct PYPLM {
DiscountResampler(const PYPLM& m) : m_(m) {}
const PYPLM& m_;
double operator()(const double& proposed_discount) const {
- return m_.log_likelihood(proposed_discount, m_.alpha);
+ return m_.log_likelihood(proposed_discount, m_.strength);
}
};
struct AlphaResampler {
AlphaResampler(const PYPLM& m) : m_(m) {}
const PYPLM& m_;
- double operator()(const double& proposed_alpha) const {
- return m_.log_likelihood(m_.d, proposed_alpha);
+ double operator()(const double& proposed_strength) const {
+ return m_.log_likelihood(m_.d, proposed_strength);
}
};
@@ -141,25 +141,25 @@ template <unsigned N> struct PYPLM {
DiscountResampler dr(*this);
AlphaResampler ar(*this);
for (int iter = 0; iter < nloop; ++iter) {
- alpha = slice_sampler1d(ar, alpha, *rng, 0.0,
+ strength = slice_sampler1d(ar, strength, *rng, 0.0,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
d = slice_sampler1d(dr, d, *rng, std::numeric_limits<double>::min(),
1.0, 0.0, niterations, 100*niterations);
}
- alpha = slice_sampler1d(ar, alpha, *rng, 0.0,
+ strength = slice_sampler1d(ar, strength, *rng, 0.0,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
typename unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > >::iterator it;
- cerr << "PYPLM<" << N << ">(d=" << d << ",a=" << alpha << ") = " << log_likelihood(d, alpha) << endl;
+ cerr << "PYPLM<" << N << ">(d=" << d << ",a=" << strength << ") = " << log_likelihood(d, strength) << endl;
for (it = p.begin(); it != p.end(); ++it) {
it->second.set_discount(d);
- it->second.set_alpha(alpha);
+ it->second.set_strength(strength);
}
backoff.resample_hyperparameters(rng, nloop, niterations);
}
PYPLM<N-1> backoff;
double discount_a, discount_b, strength_s, strength_r;
- double d, alpha;
+ double d, strength;
mutable vector<WordID> lookup; // thread-local
unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > > p;
};
diff --git a/phrasinator/gibbs_train_plm.cc b/phrasinator/gibbs_train_plm.cc
index 54861dcb..3b99e1b6 100644
--- a/phrasinator/gibbs_train_plm.cc
+++ b/phrasinator/gibbs_train_plm.cc
@@ -252,7 +252,7 @@ struct UniphraseLM {
void ResampleHyperparameters(MT19937* rng) {
phrases_.resample_hyperparameters(rng);
gen_.resample_hyperparameters(rng);
- cerr << " d=" << phrases_.discount() << ",a=" << phrases_.alpha();
+ cerr << " d=" << phrases_.discount() << ",s=" << phrases_.strength();
}
CCRP<vector<int> > phrases_;
diff --git a/utils/ccrp.h b/utils/ccrp.h
index c883c027..5f9db7a6 100644
--- a/utils/ccrp.h
+++ b/utils/ccrp.h
@@ -18,27 +18,27 @@
template <typename Dish, typename DishHash = boost::hash<Dish> >
class CCRP {
public:
- CCRP(double disc, double alpha) :
+ CCRP(double disc, double strength) :
num_tables_(),
num_customers_(),
discount_(disc),
- alpha_(alpha),
- discount_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
+ strength_(strength),
+ discount_prior_strength_(std::numeric_limits<double>::quiet_NaN()),
discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
- alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
- alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {
+ strength_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+ strength_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {
check_hyperparameters();
}
- CCRP(double d_alpha, double d_beta, double c_shape, double c_rate, double d = 0.9, double c = 1.0) :
+ CCRP(double d_strength, double d_beta, double c_shape, double c_rate, double d = 0.9, double c = 1.0) :
num_tables_(),
num_customers_(),
discount_(d),
- alpha_(c),
- discount_prior_alpha_(d_alpha),
+ strength_(c),
+ discount_prior_strength_(d_strength),
discount_prior_beta_(d_beta),
- alpha_prior_shape_(c_shape),
- alpha_prior_rate_(c_rate) {
+ strength_prior_shape_(c_shape),
+ strength_prior_rate_(c_rate) {
check_hyperparameters();
}
@@ -47,23 +47,23 @@ class CCRP {
std::cerr << "Bad discount: " << discount_ << std::endl;
abort();
}
- if (alpha_ <= -discount_) {
- std::cerr << "Bad strength: " << alpha_ << " (discount=" << discount_ << ")" << std::endl;
+ if (strength_ <= -discount_) {
+ std::cerr << "Bad strength: " << strength_ << " (discount=" << discount_ << ")" << std::endl;
abort();
}
}
double discount() const { return discount_; }
- double alpha() const { return alpha_; }
+ double strength() const { return strength_; }
void set_discount(double d) { discount_ = d; check_hyperparameters(); }
- void set_alpha(double a) { alpha_ = a; check_hyperparameters(); }
+ void set_strength(double a) { strength_ = a; check_hyperparameters(); }
bool has_discount_prior() const {
- return !std::isnan(discount_prior_alpha_);
+ return !std::isnan(discount_prior_strength_);
}
- bool has_alpha_prior() const {
- return !std::isnan(alpha_prior_shape_);
+ bool has_strength_prior() const {
+ return !std::isnan(strength_prior_shape_);
}
void clear() {
@@ -97,7 +97,7 @@ class CCRP {
DishLocations& loc = dish_locs_[dish];
bool share_table = false;
if (loc.total_dish_count_) {
- const double p_empty = (alpha_ + num_tables_ * discount_) * p0;
+ const double p_empty = (strength_ + num_tables_ * discount_) * p0;
const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
share_table = rng->SelectSample(p_empty, p_share);
}
@@ -131,7 +131,7 @@ class CCRP {
DishLocations& loc = dish_locs_[dish];
bool share_table = false;
if (loc.total_dish_count_) {
- const T p_empty = T(alpha_ + num_tables_ * discount_) * p0;
+ const T p_empty = T(strength_ + num_tables_ * discount_) * p0;
const T p_share = T(loc.total_dish_count_ - loc.table_counts_.size() * discount_);
share_table = rng->SelectSample(p_empty, p_share);
}
@@ -198,47 +198,47 @@ class CCRP {
double prob(const Dish& dish, const double& p0) const {
const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
- const double r = num_tables_ * discount_ + alpha_;
+ const double r = num_tables_ * discount_ + strength_;
if (it == dish_locs_.end()) {
- return r * p0 / (num_customers_ + alpha_);
+ return r * p0 / (num_customers_ + strength_);
} else {
return (it->second.total_dish_count_ - discount_ * it->second.table_counts_.size() + r * p0) /
- (num_customers_ + alpha_);
+ (num_customers_ + strength_);
}
}
template <typename T>
T probT(const Dish& dish, const T& p0) const {
const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
- const T r = T(num_tables_ * discount_ + alpha_);
+ const T r = T(num_tables_ * discount_ + strength_);
if (it == dish_locs_.end()) {
- return r * p0 / T(num_customers_ + alpha_);
+ return r * p0 / T(num_customers_ + strength_);
} else {
return (T(it->second.total_dish_count_ - discount_ * it->second.table_counts_.size()) + r * p0) /
- T(num_customers_ + alpha_);
+ T(num_customers_ + strength_);
}
}
double log_crp_prob() const {
- return log_crp_prob(discount_, alpha_);
+ return log_crp_prob(discount_, strength_);
}
// taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
// does not include P_0's
- double log_crp_prob(const double& discount, const double& alpha) const {
+ double log_crp_prob(const double& discount, const double& strength) const {
double lp = 0.0;
if (has_discount_prior())
- lp = Md::log_beta_density(discount, discount_prior_alpha_, discount_prior_beta_);
- if (has_alpha_prior())
- lp += Md::log_gamma_density(alpha + discount, alpha_prior_shape_, alpha_prior_rate_);
+ lp = Md::log_beta_density(discount, discount_prior_strength_, discount_prior_beta_);
+ if (has_strength_prior())
+ lp += Md::log_gamma_density(strength + discount, strength_prior_shape_, strength_prior_rate_);
assert(lp <= 0.0);
if (num_customers_) {
if (discount > 0.0) {
const double r = lgamma(1.0 - discount);
- if (alpha)
- lp += lgamma(alpha) - lgamma(alpha / discount);
- lp += - lgamma(alpha + num_customers_)
- + num_tables_ * log(discount) + lgamma(alpha / discount + num_tables_);
+ if (strength)
+ lp += lgamma(strength) - lgamma(strength / discount);
+ lp += - lgamma(strength + num_customers_)
+ + num_tables_ * log(discount) + lgamma(strength / discount + num_tables_);
assert(std::isfinite(lp));
for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
it != dish_locs_.end(); ++it) {
@@ -247,8 +247,16 @@ class CCRP {
lp += lgamma(*ti - discount) - r;
}
}
+ } else if (!discount) { // discount == 0.0
+ lp += lgamma(strength) + num_tables_ * log(strength) - lgamma(strength + num_tables_);
+ assert(std::isfinite(lp));
+ for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
+ it != dish_locs_.end(); ++it) {
+ const DishLocations& cur = it->second;
+ lp += lgamma(cur.table_counts_.size());
+ }
} else {
- assert(!"not implemented yet");
+ assert(!"discount less than 0 detected!");
}
}
assert(std::isfinite(lp));
@@ -256,22 +264,22 @@ class CCRP {
}
void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
- assert(has_discount_prior() || has_alpha_prior());
+ assert(has_discount_prior() || has_strength_prior());
DiscountResampler dr(*this);
StrengthResampler sr(*this);
for (int iter = 0; iter < nloop; ++iter) {
- if (has_alpha_prior()) {
- alpha_ = slice_sampler1d(sr, alpha_, *rng, -discount_,
+ if (has_strength_prior()) {
+ strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
if (has_discount_prior()) {
double min_discount = std::numeric_limits<double>::min();
- if (alpha_ < 0.0) min_discount = -alpha_;
+ if (strength_ < 0.0) min_discount = -strength_;
discount_ = slice_sampler1d(dr, discount_, *rng, min_discount,
1.0, 0.0, niterations, 100*niterations);
}
}
- alpha_ = slice_sampler1d(sr, alpha_, *rng, -discount_,
+ strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
@@ -279,15 +287,15 @@ class CCRP {
DiscountResampler(const CCRP& crp) : crp_(crp) {}
const CCRP& crp_;
double operator()(const double& proposed_discount) const {
- return crp_.log_crp_prob(proposed_discount, crp_.alpha_);
+ return crp_.log_crp_prob(proposed_discount, crp_.strength_);
}
};
struct StrengthResampler {
StrengthResampler(const CCRP& crp) : crp_(crp) {}
const CCRP& crp_;
- double operator()(const double& proposed_alpha) const {
- return crp_.log_crp_prob(crp_.discount_, proposed_alpha);
+ double operator()(const double& proposed_strength) const {
+ return crp_.log_crp_prob(crp_.discount_, proposed_strength);
}
};
@@ -299,7 +307,7 @@ class CCRP {
};
void Print(std::ostream* out) const {
- std::cerr << "PYP(d=" << discount_ << ",c=" << alpha_ << ") customers=" << num_customers_ << std::endl;
+ std::cerr << "PYP(d=" << discount_ << ",c=" << strength_ << ") customers=" << num_customers_ << std::endl;
for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
it != dish_locs_.end(); ++it) {
(*out) << it->first << " (" << it->second.total_dish_count_ << " on " << it->second.table_counts_.size() << " tables): ";
@@ -324,15 +332,15 @@ class CCRP {
std::tr1::unordered_map<Dish, DishLocations, DishHash> dish_locs_;
double discount_;
- double alpha_;
+ double strength_;
// optional beta prior on discount_ (NaN if no prior)
- double discount_prior_alpha_;
+ double discount_prior_strength_;
double discount_prior_beta_;
- // optional gamma prior on alpha_ (NaN if no prior)
- double alpha_prior_shape_;
- double alpha_prior_rate_;
+ // optional gamma prior on strength_ (NaN if no prior)
+ double strength_prior_shape_;
+ double strength_prior_rate_;
};
template <typename T,typename H>
diff --git a/utils/mfcr.h b/utils/mfcr.h
index df988f51..aeaf599d 100644
--- a/utils/mfcr.h
+++ b/utils/mfcr.h
@@ -39,37 +39,37 @@ template <typename Dish, typename DishHash = boost::hash<Dish> >
class MFCR {
public:
- MFCR(unsigned num_floors, double d, double alpha) :
+ MFCR(unsigned num_floors, double d, double strength) :
num_floors_(num_floors),
num_tables_(),
num_customers_(),
discount_(d),
- alpha_(alpha),
- discount_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
+ strength_(strength),
+ discount_prior_strength_(std::numeric_limits<double>::quiet_NaN()),
discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
- alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
- alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+ strength_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+ strength_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
- MFCR(unsigned num_floors, double discount_alpha, double discount_beta, double alpha_shape, double alpha_rate, double d = 0.9, double alpha = 10.0) :
+ MFCR(unsigned num_floors, double discount_strength, double discount_beta, double strength_shape, double strength_rate, double d = 0.9, double strength = 10.0) :
num_floors_(num_floors),
num_tables_(),
num_customers_(),
discount_(d),
- alpha_(alpha),
- discount_prior_alpha_(discount_alpha),
+ strength_(strength),
+ discount_prior_strength_(discount_strength),
discount_prior_beta_(discount_beta),
- alpha_prior_shape_(alpha_shape),
- alpha_prior_rate_(alpha_rate) {}
+ strength_prior_shape_(strength_shape),
+ strength_prior_rate_(strength_rate) {}
double discount() const { return discount_; }
- double alpha() const { return alpha_; }
+ double strength() const { return strength_; }
bool has_discount_prior() const {
- return !std::isnan(discount_prior_alpha_);
+ return !std::isnan(discount_prior_strength_);
}
- bool has_alpha_prior() const {
- return !std::isnan(alpha_prior_shape_);
+ bool has_strength_prior() const {
+ return !std::isnan(strength_prior_shape_);
}
void clear() {
@@ -122,7 +122,7 @@ class MFCR {
int floor = -1;
bool share_table = false;
if (loc.total_dish_count_) {
- const double p_empty = (alpha_ + num_tables_ * discount_) * marg_p0;
+ const double p_empty = (strength_ + num_tables_ * discount_) * marg_p0;
const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
share_table = rng->SelectSample(p_empty, p_share);
}
@@ -206,44 +206,53 @@ class MFCR {
const double marg_p0 = std::inner_product(p0s.begin(), p0s.end(), lambdas.begin(), 0.0);
assert(marg_p0 <= 1.0);
const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
- const double r = num_tables_ * discount_ + alpha_;
+ const double r = num_tables_ * discount_ + strength_;
if (it == dish_locs_.end()) {
- return r * marg_p0 / (num_customers_ + alpha_);
+ return r * marg_p0 / (num_customers_ + strength_);
} else {
return (it->second.total_dish_count_ - discount_ * it->second.table_counts_.size() + r * marg_p0) /
- (num_customers_ + alpha_);
+ (num_customers_ + strength_);
}
}
double log_crp_prob() const {
- return log_crp_prob(discount_, alpha_);
+ return log_crp_prob(discount_, strength_);
}
// taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
// does not include draws from G_w's
- double log_crp_prob(const double& d, const double& alpha) const {
+ double log_crp_prob(const double& discount, const double& strength) const {
double lp = 0.0;
if (has_discount_prior())
- lp = Md::log_beta_density(d, discount_prior_alpha_, discount_prior_beta_);
- if (has_alpha_prior())
- lp += Md::log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
+ lp = Md::log_beta_density(discount, discount_prior_strength_, discount_prior_beta_);
+ if (has_strength_prior())
+ lp += Md::log_gamma_density(strength + discount, strength_prior_shape_, strength_prior_rate_);
assert(lp <= 0.0);
if (num_customers_) {
- if (d > 0.0) {
- const double r = lgamma(1.0 - d);
- lp += lgamma(alpha) - lgamma(alpha + num_customers_)
- + num_tables_ * log(d) + lgamma(alpha / d + num_tables_)
- - lgamma(alpha / d);
+ if (discount > 0.0) {
+ const double r = lgamma(1.0 - discount);
+ if (strength)
+ lp += lgamma(strength) - lgamma(strength / discount);
+ lp += - lgamma(strength + num_customers_)
+ + num_tables_ * log(discount) + lgamma(strength / discount + num_tables_);
assert(std::isfinite(lp));
for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
it != dish_locs_.end(); ++it) {
const DishLocations& cur = it->second;
for (std::list<TableCount>::const_iterator ti = cur.table_counts_.begin(); ti != cur.table_counts_.end(); ++ti) {
- lp += lgamma(ti->count - d) - r;
+ lp += lgamma(ti->count - discount) - r;
}
}
+ } else if (!discount) { // discount == 0.0
+ lp += lgamma(strength) + num_tables_ * log(strength) - lgamma(strength + num_tables_);
+ assert(std::isfinite(lp));
+ for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
+ it != dish_locs_.end(); ++it) {
+ const DishLocations& cur = it->second;
+ lp += lgamma(cur.table_counts_.size());
+ }
} else {
- assert(!"not implemented yet");
+ assert(!"discount less than 0 detected!");
}
}
assert(std::isfinite(lp));
@@ -251,20 +260,22 @@ class MFCR {
}
void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
- assert(has_discount_prior() || has_alpha_prior());
+ assert(has_discount_prior() || has_strength_prior());
DiscountResampler dr(*this);
- ConcentrationResampler cr(*this);
+ StrengthResampler sr(*this);
for (int iter = 0; iter < nloop; ++iter) {
- if (has_alpha_prior()) {
- alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
+ if (has_strength_prior()) {
+ strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
if (has_discount_prior()) {
- discount_ = slice_sampler1d(dr, discount_, *rng, std::numeric_limits<double>::min(),
+ double min_discount = std::numeric_limits<double>::min();
+ if (strength_ < 0.0) min_discount = -strength_;
+ discount_ = slice_sampler1d(dr, discount_, *rng, min_discount,
1.0, 0.0, niterations, 100*niterations);
}
}
- alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
+ strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
@@ -272,15 +283,15 @@ class MFCR {
DiscountResampler(const MFCR& crp) : crp_(crp) {}
const MFCR& crp_;
double operator()(const double& proposed_d) const {
- return crp_.log_crp_prob(proposed_d, crp_.alpha_);
+ return crp_.log_crp_prob(proposed_d, crp_.strength_);
}
};
- struct ConcentrationResampler {
- ConcentrationResampler(const MFCR& crp) : crp_(crp) {}
+ struct StrengthResampler {
+ StrengthResampler(const MFCR& crp) : crp_(crp) {}
const MFCR& crp_;
- double operator()(const double& proposediscount_alpha) const {
- return crp_.log_crp_prob(crp_.discount_, proposediscount_alpha);
+ double operator()(const double& proposediscount_strength) const {
+ return crp_.log_crp_prob(crp_.discount_, proposediscount_strength);
}
};
@@ -292,7 +303,7 @@ class MFCR {
};
void Print(std::ostream* out) const {
- (*out) << "MFCR(d=" << discount_ << ",alpha=" << alpha_ << ") customers=" << num_customers_ << std::endl;
+ (*out) << "MFCR(d=" << discount_ << ",strength=" << strength_ << ") customers=" << num_customers_ << std::endl;
for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
it != dish_locs_.end(); ++it) {
(*out) << it->first << " (" << it->second.total_dish_count_ << " on " << it->second.table_counts_.size() << " tables): ";
@@ -318,15 +329,15 @@ class MFCR {
std::tr1::unordered_map<Dish, DishLocations, DishHash> dish_locs_;
double discount_;
- double alpha_;
+ double strength_;
// optional beta prior on discount_ (NaN if no prior)
- double discount_prior_alpha_;
+ double discount_prior_strength_;
double discount_prior_beta_;
- // optional gamma prior on alpha_ (NaN if no prior)
- double alpha_prior_shape_;
- double alpha_prior_rate_;
+ // optional gamma prior on strength_ (NaN if no prior)
+ double strength_prior_shape_;
+ double strength_prior_rate_;
};
template <typename T,typename H>
--
cgit v1.2.3
From 4c007d48d5829233d0ae3c3c8b48f8c25631bf81 Mon Sep 17 00:00:00 2001
From: Chris Dyer <cdyer@cs.cmu.edu>
Date: Mon, 5 Mar 2012 16:06:45 -0500
Subject: use template parameter inference to figure out what type to use for
probability computations, templatatize number of floors in MFCR rather than
compile-time set
---
gi/pf/align-lexonly-pyp.cc | 20 +++++++-------
gi/pf/conditional_pseg.h | 22 +++++++--------
gi/pf/learn_cfg.cc | 8 +++---
utils/ccrp.h | 48 ++------------------------------
utils/mfcr.h | 68 ++++++++++++++++++++++++----------------------
utils/mfcr_test.cc | 10 +++----
6 files changed, 68 insertions(+), 108 deletions(-)
(limited to 'utils/mfcr.h')
diff --git a/gi/pf/align-lexonly-pyp.cc b/gi/pf/align-lexonly-pyp.cc
index 87f7f6b5..ac0590e0 100644
--- a/gi/pf/align-lexonly-pyp.cc
+++ b/gi/pf/align-lexonly-pyp.cc
@@ -68,7 +68,7 @@ struct AlignedSentencePair {
struct HierarchicalWordBase {
explicit HierarchicalWordBase(const unsigned vocab_e_size) :
- base(prob_t::One()), r(1,1,1,25,25), u0(-log(vocab_e_size)), l(1,1.0), v(1, 0.0) {}
+ base(prob_t::One()), r(1,1,1,1), u0(-log(vocab_e_size)), l(1,prob_t::One()), v(1, prob_t::Zero()) {}
void ResampleHyperparameters(MT19937* rng) {
r.resample_hyperparameters(rng);
@@ -80,14 +80,14 @@ struct HierarchicalWordBase {
// return p0 of rule.e_
prob_t operator()(const TRule& rule) const {
- v[0] = exp(logp0(rule.e_));
- return prob_t(r.prob(rule.e_, v, l));
+ v[0].logeq(logp0(rule.e_));
+ return r.prob(rule.e_, v.begin(), l.begin());
}
void Increment(const TRule& rule) {
- v[0] = exp(logp0(rule.e_));
- if (r.increment(rule.e_, v, l, &*prng).count) {
- base *= prob_t(v[0] * l[0]);
+ v[0].logeq(logp0(rule.e_));
+ if (r.increment(rule.e_, v.begin(), l.begin(), &*prng).count) {
+ base *= v[0] * l[0];
}
}
@@ -105,15 +105,15 @@ struct HierarchicalWordBase {
void Summary() const {
cerr << "NUMBER OF CUSTOMERS: " << r.num_customers() << " (d=" << r.discount() << ",s=" << r.strength() << ')' << endl;
- for (MFCR<vector<WordID> >::const_iterator it = r.begin(); it != r.end(); ++it)
+ for (MFCR<1,vector<WordID> >::const_iterator it = r.begin(); it != r.end(); ++it)
cerr << " " << it->second.total_dish_count_ << " (on " << it->second.table_counts_.size() << " tables)" << TD::GetString(it->first) << endl;
}
prob_t base;
- MFCR<vector<WordID> > r;
+ MFCR<1,vector<WordID> > r;
const double u0;
- const vector<double> l;
- mutable vector<double> v;
+ const vector<prob_t> l;
+ mutable vector<prob_t> v;
};
struct BasicLexicalAlignment {
diff --git a/gi/pf/conditional_pseg.h b/gi/pf/conditional_pseg.h
index 86403d8d..ef73e332 100644
--- a/gi/pf/conditional_pseg.h
+++ b/gi/pf/conditional_pseg.h
@@ -17,13 +17,13 @@
template <typename ConditionalBaseMeasure>
struct MConditionalTranslationModel {
explicit MConditionalTranslationModel(ConditionalBaseMeasure& rcp0) :
- rp0(rcp0), lambdas(1, 1.0), p0s(1) {}
+ rp0(rcp0), lambdas(1, prob_t::One()), p0s(1) {}
void Summary() const {
std::cerr << "Number of conditioning contexts: " << r.size() << std::endl;
for (RuleModelHash::const_iterator it = r.begin(); it != r.end(); ++it) {
std::cerr << TD::GetString(it->first) << " \t(d=" << it->second.discount() << ",s=" << it->second.strength() << ") --------------------------" << std::endl;
- for (MFCR<TRule>::const_iterator i2 = it->second.begin(); i2 != it->second.end(); ++i2)
+ for (MFCR<1,TRule>::const_iterator i2 = it->second.begin(); i2 != it->second.end(); ++i2)
std::cerr << " " << -1 << '\t' << i2->first << std::endl;
}
}
@@ -46,10 +46,10 @@ struct MConditionalTranslationModel {
int IncrementRule(const TRule& rule, MT19937* rng) {
RuleModelHash::iterator it = r.find(rule.f_);
if (it == r.end()) {
- it = r.insert(make_pair(rule.f_, MFCR<TRule>(1, 1.0, 1.0, 1.0, 1.0, 1e-9, 4.0))).first;
+ it = r.insert(make_pair(rule.f_, MFCR<1,TRule>(1.0, 1.0, 1.0, 1.0, 1e-9, 4.0))).first;
}
- p0s[0] = rp0(rule).as_float();
- TableCount delta = it->second.increment(rule, p0s, lambdas, rng);
+ p0s[0] = rp0(rule);
+ TableCount delta = it->second.increment(rule, p0s.begin(), lambdas.begin(), rng);
return delta.count;
}
@@ -57,10 +57,10 @@ struct MConditionalTranslationModel {
prob_t p;
RuleModelHash::const_iterator it = r.find(rule.f_);
if (it == r.end()) {
- p.logeq(log(rp0(rule)));
+ p = rp0(rule);
} else {
- p0s[0] = rp0(rule).as_float();
- p = prob_t(it->second.prob(rule, p0s, lambdas));
+ p0s[0] = rp0(rule);
+ p = it->second.prob(rule, p0s.begin(), lambdas.begin());
}
return p;
}
@@ -80,11 +80,11 @@ struct MConditionalTranslationModel {
const ConditionalBaseMeasure& rp0;
typedef std::tr1::unordered_map<std::vector<WordID>,
- MFCR<TRule>,
+ MFCR<1, TRule>,
boost::hash<std::vector<WordID> > > RuleModelHash;
RuleModelHash r;
- std::vector<double> lambdas;
- mutable std::vector<double> p0s;
+ std::vector<prob_t> lambdas;
+ mutable std::vector<prob_t> p0s;
};
template <typename ConditionalBaseMeasure>
diff --git a/gi/pf/learn_cfg.cc b/gi/pf/learn_cfg.cc
index bf157828..ed1772bf 100644
--- a/gi/pf/learn_cfg.cc
+++ b/gi/pf/learn_cfg.cc
@@ -127,20 +127,20 @@ struct HieroLMModel {
nts(num_nts, CCRP<TRule>(1,1,1,1)) {}
prob_t Prob(const TRule& r) const {
- return nts[nt_id_to_index[-r.lhs_]].probT<prob_t>(r, p0(r));
+ return nts[nt_id_to_index[-r.lhs_]].prob(r, p0(r));
}
inline prob_t p0(const TRule& r) const {
if (kHIERARCHICAL_PRIOR)
- return q0.probT<prob_t>(r, base(r));
+ return q0.prob(r, base(r));
else
return base(r);
}
int Increment(const TRule& r, MT19937* rng) {
- const int delta = nts[nt_id_to_index[-r.lhs_]].incrementT<prob_t>(r, p0(r), rng);
+ const int delta = nts[nt_id_to_index[-r.lhs_]].increment(r, p0(r), rng);
if (kHIERARCHICAL_PRIOR && delta)
- q0.incrementT<prob_t>(r, base(r), rng);
+ q0.increment(r, base(r), rng);
return delta;
// return x.increment(r);
}
diff --git a/utils/ccrp.h b/utils/ccrp.h
index 5f9db7a6..e24130ac 100644
--- a/utils/ccrp.h
+++ b/utils/ccrp.h
@@ -92,42 +92,9 @@ class CCRP {
return it->total_dish_count_;
}
- // returns +1 or 0 indicating whether a new table was opened
- int increment(const Dish& dish, const double& p0, MT19937* rng) {
- DishLocations& loc = dish_locs_[dish];
- bool share_table = false;
- if (loc.total_dish_count_) {
- const double p_empty = (strength_ + num_tables_ * discount_) * p0;
- const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
- share_table = rng->SelectSample(p_empty, p_share);
- }
- if (share_table) {
- double r = rng->next() * (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
- for (typename std::list<unsigned>::iterator ti = loc.table_counts_.begin();
- ti != loc.table_counts_.end(); ++ti) {
- r -= (*ti - discount_);
- if (r <= 0.0) {
- ++(*ti);
- break;
- }
- }
- if (r > 0.0) {
- std::cerr << "Serious error: r=" << r << std::endl;
- Print(&std::cerr);
- assert(r <= 0.0);
- }
- } else {
- loc.table_counts_.push_back(1u);
- ++num_tables_;
- }
- ++loc.total_dish_count_;
- ++num_customers_;
- return (share_table ? 0 : 1);
- }
-
// returns +1 or 0 indicating whether a new table was opened
template <typename T>
- int incrementT(const Dish& dish, const T& p0, MT19937* rng) {
+ int increment(const Dish& dish, const T& p0, MT19937* rng) {
DishLocations& loc = dish_locs_[dish];
bool share_table = false;
if (loc.total_dish_count_) {
@@ -196,19 +163,8 @@ class CCRP {
}
}
- double prob(const Dish& dish, const double& p0) const {
- const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
- const double r = num_tables_ * discount_ + strength_;
- if (it == dish_locs_.end()) {
- return r * p0 / (num_customers_ + strength_);
- } else {
- return (it->second.total_dish_count_ - discount_ * it->second.table_counts_.size() + r * p0) /
- (num_customers_ + strength_);
- }
- }
-
template <typename T>
- T probT(const Dish& dish, const T& p0) const {
+ T prob(const Dish& dish, const T& p0) const {
const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
const T r = T(num_tables_ * discount_ + strength_);
if (it == dish_locs_.end()) {
diff --git a/utils/mfcr.h b/utils/mfcr.h
index aeaf599d..6cc0ebf1 100644
--- a/utils/mfcr.h
+++ b/utils/mfcr.h
@@ -8,6 +8,7 @@
#include <list>
#include <iostream>
#include <vector>
+#include <iterator>
#include <tr1/unordered_map>
#include <boost/functional/hash.hpp>
#include "sampler.h"
@@ -35,12 +36,11 @@ std::ostream& operator<<(std::ostream& o, const TableCount& tc) {
// referenced therein.
// http://www.aclweb.org/anthology/P/P09/P09-2085.pdf
//
-template <typename Dish, typename DishHash = boost::hash<Dish> >
+template <unsigned Floors, typename Dish, typename DishHash = boost::hash<Dish> >
class MFCR {
public:
- MFCR(unsigned num_floors, double d, double strength) :
- num_floors_(num_floors),
+ MFCR(double d, double strength) :
num_tables_(),
num_customers_(),
discount_(d),
@@ -50,8 +50,7 @@ class MFCR {
strength_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
strength_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
- MFCR(unsigned num_floors, double discount_strength, double discount_beta, double strength_shape, double strength_rate, double d = 0.9, double strength = 10.0) :
- num_floors_(num_floors),
+ MFCR(double discount_strength, double discount_beta, double strength_shape, double strength_rate, double d = 0.9, double strength = 10.0) :
num_tables_(),
num_customers_(),
discount_(d),
@@ -111,22 +110,22 @@ class MFCR {
}
// returns (delta, floor) indicating whether a new table (delta) was opened and on which floor
- TableCount increment(const Dish& dish, const std::vector<double>& p0s, const std::vector<double>& lambdas, MT19937* rng) {
- assert(p0s.size() == num_floors_);
- assert(lambdas.size() == num_floors_);
-
+ template <class InputIterator, class InputIterator2>
+ TableCount increment(const Dish& dish, InputIterator p0s, InputIterator2 lambdas, MT19937* rng) {
DishLocations& loc = dish_locs_[dish];
// marg_p0 = marginal probability of opening a new table on any floor with label dish
- const double marg_p0 = std::inner_product(p0s.begin(), p0s.end(), lambdas.begin(), 0.0);
- assert(marg_p0 <= 1.0);
+ typedef typename std::iterator_traits<InputIterator>::value_type F;
+ const F marg_p0 = std::inner_product(p0s, p0s + Floors, lambdas, F(0.0));
+ assert(marg_p0 <= F(1.0001));
int floor = -1;
bool share_table = false;
if (loc.total_dish_count_) {
- const double p_empty = (strength_ + num_tables_ * discount_) * marg_p0;
- const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
+ const F p_empty = F(strength_ + num_tables_ * discount_) * marg_p0;
+ const F p_share = F(loc.total_dish_count_ - loc.table_counts_.size() * discount_);
share_table = rng->SelectSample(p_empty, p_share);
}
if (share_table) {
+ // this can be done with doubles since P0 (which may be tiny) is not involved
double r = rng->next() * (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
for (typename std::list<TableCount>::iterator ti = loc.table_counts_.begin();
ti != loc.table_counts_.end(); ++ti) {
@@ -143,12 +142,18 @@ class MFCR {
assert(r <= 0.0);
}
} else { // sit at currently empty table -- must sample what floor
- double r = rng->next() * marg_p0;
- for (unsigned i = 0; i < p0s.size(); ++i) {
- r -= p0s[i] * lambdas[i];
- if (r <= 0.0) {
- floor = i;
- break;
+ if (Floors == 1) {
+ floor = 0;
+ } else {
+ F r = F(rng->next()) * marg_p0;
+ for (unsigned i = 0; i < Floors; ++i) {
+ r -= (*p0s) * (*lambdas);
+ ++p0s;
+ ++lambdas;
+ if (r <= F(0.0)) {
+ floor = i;
+ break;
+ }
}
}
assert(floor >= 0);
@@ -200,18 +205,18 @@ class MFCR {
return TableCount(delta, floor);
}
- double prob(const Dish& dish, const std::vector<double>& p0s, const std::vector<double>& lambdas) const {
- assert(p0s.size() == num_floors_);
- assert(lambdas.size() == num_floors_);
- const double marg_p0 = std::inner_product(p0s.begin(), p0s.end(), lambdas.begin(), 0.0);
- assert(marg_p0 <= 1.0);
+ template <class InputIterator, class InputIterator2>
+ typename std::iterator_traits<InputIterator>::value_type prob(const Dish& dish, InputIterator p0s, InputIterator2 lambdas) const {
+ typedef typename std::iterator_traits<InputIterator>::value_type F;
+ const F marg_p0 = std::inner_product(p0s, p0s + Floors, lambdas, F(0.0));
+ assert(marg_p0 <= F(1.0001));
const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
- const double r = num_tables_ * discount_ + strength_;
+ const F r = F(num_tables_ * discount_ + strength_);
if (it == dish_locs_.end()) {
- return r * marg_p0 / (num_customers_ + strength_);
+ return r * marg_p0 / F(num_customers_ + strength_);
} else {
- return (it->second.total_dish_count_ - discount_ * it->second.table_counts_.size() + r * marg_p0) /
- (num_customers_ + strength_);
+ return (F(it->second.total_dish_count_ - discount_ * it->second.table_counts_.size()) + F(r * marg_p0)) /
+ F(num_customers_ + strength_);
}
}
@@ -303,7 +308,7 @@ class MFCR {
};
void Print(std::ostream* out) const {
- (*out) << "MFCR(d=" << discount_ << ",strength=" << strength_ << ") customers=" << num_customers_ << std::endl;
+ (*out) << "MFCR<" << Floors << ">(d=" << discount_ << ",strength=" << strength_ << ") customers=" << num_customers_ << std::endl;
for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
it != dish_locs_.end(); ++it) {
(*out) << it->first << " (" << it->second.total_dish_count_ << " on " << it->second.table_counts_.size() << " tables): ";
@@ -323,7 +328,6 @@ class MFCR {
return dish_locs_.end();
}
- unsigned num_floors_;
unsigned num_tables_;
unsigned num_customers_;
std::tr1::unordered_map<Dish, DishLocations, DishHash> dish_locs_;
@@ -340,8 +344,8 @@ class MFCR {
double strength_prior_rate_;
};
-template <typename T,typename H>
-std::ostream& operator<<(std::ostream& o, const MFCR<T,H>& c) {
+template <unsigned N,typename T,typename H>
+std::ostream& operator<<(std::ostream& o, const MFCR<N,T,H>& c) {
c.Print(&o);
return o;
}
diff --git a/utils/mfcr_test.cc b/utils/mfcr_test.cc
index 7c45a37c..cc886335 100644
--- a/utils/mfcr_test.cc
+++ b/utils/mfcr_test.cc
@@ -9,7 +9,7 @@
using namespace std;
void test_exch(MT19937* rng) {
- MFCR<int> crp(2, 0.5, 3.0);
+ MFCR<2, int> crp(0.5, 3.0);
vector<double> lambdas(2);
vector<double> p0s(2);
lambdas[0] = 0.2;
@@ -22,23 +22,23 @@ void test_exch(MT19937* rng) {
double xt = 0;
int cust = 10;
vector<int> hist(cust + 1, 0), hist2(cust + 1, 0);
- for (int i = 0; i < cust; ++i) { crp.increment(1, p0s, lambdas, rng); }
+ for (int i = 0; i < cust; ++i) { crp.increment(1, p0s.begin(), lambdas.begin(), rng); }
const int samples = 100000;
const bool simulate = true;
for (int k = 0; k < samples; ++k) {
if (!simulate) {
crp.clear();
- for (int i = 0; i < cust; ++i) { crp.increment(1, p0s, lambdas, rng); }
+ for (int i = 0; i < cust; ++i) { crp.increment(1, p0s.begin(), lambdas.begin(), rng); }
} else {
int da = rng->next() * cust;
bool a = rng->next() < 0.45;
if (a) {
- for (int i = 0; i < da; ++i) { crp.increment(1, p0s, lambdas, rng); }
+ for (int i = 0; i < da; ++i) { crp.increment(1, p0s.begin(), lambdas.begin(), rng); }
for (int i = 0; i < da; ++i) { crp.decrement(1, rng); }
xt += 1.0;
} else {
for (int i = 0; i < da; ++i) { crp.decrement(1, rng); }
- for (int i = 0; i < da; ++i) { crp.increment(1, p0s, lambdas, rng); }
+ for (int i = 0; i < da; ++i) { crp.increment(1, p0s.begin(), lambdas.begin(), rng); }
}
}
int c = crp.num_tables(1);
--
cgit v1.2.3
From 7b3936660fb777b455079c63c23aec00f60f98ea Mon Sep 17 00:00:00 2001
From: Chris Dyer <cdyer@cs.cmu.edu>
Date: Mon, 5 Mar 2012 21:36:07 -0500
Subject: tie hyperparameters for translation distributions; support theta < 0
for PYPLM
---
gi/pf/align-lexonly-pyp.cc | 13 ++++-----
gi/pf/conditional_pseg.h | 68 ++++++++++++++++++++++++++++++++++++----------
gi/pf/pyp_lm.cc | 12 ++++----
utils/ccrp.h | 4 +--
utils/mfcr.h | 19 +++++++++++--
5 files changed, 84 insertions(+), 32 deletions(-)
(limited to 'utils/mfcr.h')
diff --git a/gi/pf/align-lexonly-pyp.cc b/gi/pf/align-lexonly-pyp.cc
index ac0590e0..13a3a487 100644
--- a/gi/pf/align-lexonly-pyp.cc
+++ b/gi/pf/align-lexonly-pyp.cc
@@ -68,14 +68,14 @@ struct AlignedSentencePair {
struct HierarchicalWordBase {
explicit HierarchicalWordBase(const unsigned vocab_e_size) :
- base(prob_t::One()), r(1,1,1,1), u0(-log(vocab_e_size)), l(1,prob_t::One()), v(1, prob_t::Zero()) {}
+ base(prob_t::One()), r(1,1,1,1,0.66,50.0), u0(-log(vocab_e_size)), l(1,prob_t::One()), v(1, prob_t::Zero()) {}
void ResampleHyperparameters(MT19937* rng) {
r.resample_hyperparameters(rng);
}
inline double logp0(const vector<WordID>& s) const {
- return s.size() * u0;
+ return Md::log_poisson(s.size(), 7.5) + s.size() * u0;
}
// return p0 of rule.e_
@@ -106,7 +106,7 @@ struct HierarchicalWordBase {
void Summary() const {
cerr << "NUMBER OF CUSTOMERS: " << r.num_customers() << " (d=" << r.discount() << ",s=" << r.strength() << ')' << endl;
for (MFCR<1,vector<WordID> >::const_iterator it = r.begin(); it != r.end(); ++it)
- cerr << " " << it->second.total_dish_count_ << " (on " << it->second.table_counts_.size() << " tables)" << TD::GetString(it->first) << endl;
+ cerr << " " << it->second.total_dish_count_ << " (on " << it->second.table_counts_.size() << " tables) " << TD::GetString(it->first) << endl;
}
prob_t base;
@@ -167,10 +167,9 @@ struct BasicLexicalAlignment {
}
void ResampleHyperparemeters() {
- cerr << " LLH_prev = " << Likelihood() << flush;
tmodel.ResampleHyperparameters(&*prng);
up0.ResampleHyperparameters(&*prng);
- cerr << "\tLLH_post = " << Likelihood() << endl;
+ cerr << " (base d=" << up0.r.discount() << ",s=" << up0.r.strength() << ")\n";
}
void ResampleCorpus();
@@ -218,7 +217,7 @@ void BasicLexicalAlignment::ResampleCorpus() {
up0.Increment(r);
}
}
- cerr << " LLH = " << tmodel.Likelihood() << endl;
+ cerr << " LLH = " << Likelihood() << endl;
}
void ExtractLetters(const set<WordID>& v, vector<vector<WordID> >* l, set<WordID>* letset = NULL) {
@@ -311,7 +310,7 @@ int main(int argc, char** argv) {
for (int i = 0; i < samples; ++i) {
for (int j = 65; j < 67; ++j) Debug(corpus[j]);
cerr << i << "\t" << x.tmodel.r.size() << "\t";
- if (i % 10 == 0) x.ResampleHyperparemeters();
+ if (i % 7 == 6) x.ResampleHyperparemeters();
x.ResampleCorpus();
if (i > (samples / 5) && (i % 10 == 9)) for (int j = 0; j < corpus.size(); ++j) AddSample(&corpus[j]);
}
diff --git a/gi/pf/conditional_pseg.h b/gi/pf/conditional_pseg.h
index ef73e332..8202778b 100644
--- a/gi/pf/conditional_pseg.h
+++ b/gi/pf/conditional_pseg.h
@@ -17,21 +17,66 @@
template <typename ConditionalBaseMeasure>
struct MConditionalTranslationModel {
explicit MConditionalTranslationModel(ConditionalBaseMeasure& rcp0) :
- rp0(rcp0), lambdas(1, prob_t::One()), p0s(1) {}
+ rp0(rcp0), d(0.5), strength(1.0), lambdas(1, prob_t::One()), p0s(1) {}
void Summary() const {
std::cerr << "Number of conditioning contexts: " << r.size() << std::endl;
for (RuleModelHash::const_iterator it = r.begin(); it != r.end(); ++it) {
std::cerr << TD::GetString(it->first) << " \t(d=" << it->second.discount() << ",s=" << it->second.strength() << ") --------------------------" << std::endl;
for (MFCR<1,TRule>::const_iterator i2 = it->second.begin(); i2 != it->second.end(); ++i2)
- std::cerr << " " << -1 << '\t' << i2->first << std::endl;
+ std::cerr << " " << i2->second.total_dish_count_ << '\t' << i2->first << std::endl;
}
}
+ double log_likelihood(const double& dd, const double& aa) const {
+ if (aa <= -dd) return -std::numeric_limits<double>::infinity();
+ //double llh = Md::log_beta_density(dd, 10, 3) + Md::log_gamma_density(aa, 1, 1);
+ double llh = Md::log_beta_density(dd, 1, 1) +
+ Md::log_gamma_density(dd + aa, 1, 1);
+ typename std::tr1::unordered_map<std::vector<WordID>, MFCR<1,TRule>, boost::hash<std::vector<WordID> > >::const_iterator it;
+ for (it = r.begin(); it != r.end(); ++it)
+ llh += it->second.log_crp_prob(dd, aa);
+ return llh;
+ }
+
+ struct DiscountResampler {
+ DiscountResampler(const MConditionalTranslationModel& m) : m_(m) {}
+ const MConditionalTranslationModel& m_;
+ double operator()(const double& proposed_discount) const {
+ return m_.log_likelihood(proposed_discount, m_.strength);
+ }
+ };
+
+ struct AlphaResampler {
+ AlphaResampler(const MConditionalTranslationModel& m) : m_(m) {}
+ const MConditionalTranslationModel& m_;
+ double operator()(const double& proposed_strength) const {
+ return m_.log_likelihood(m_.d, proposed_strength);
+ }
+ };
+
void ResampleHyperparameters(MT19937* rng) {
- for (RuleModelHash::iterator it = r.begin(); it != r.end(); ++it)
- it->second.resample_hyperparameters(rng);
- }
+ const unsigned nloop = 5;
+ const unsigned niterations = 10;
+ DiscountResampler dr(*this);
+ AlphaResampler ar(*this);
+ for (int iter = 0; iter < nloop; ++iter) {
+ strength = slice_sampler1d(ar, strength, *rng, -d + std::numeric_limits<double>::min(),
+ std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+ double min_discount = std::numeric_limits<double>::min();
+ if (strength < 0.0) min_discount -= strength;
+ d = slice_sampler1d(dr, d, *rng, min_discount,
+ 1.0, 0.0, niterations, 100*niterations);
+ }
+ strength = slice_sampler1d(ar, strength, *rng, -d,
+ std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+ typename std::tr1::unordered_map<std::vector<WordID>, MFCR<1,TRule>, boost::hash<std::vector<WordID> > >::iterator it;
+ std::cerr << "MConditionalTranslationModel(d=" << d << ",s=" << strength << ") = " << log_likelihood(d, strength) << std::endl;
+ for (it = r.begin(); it != r.end(); ++it) {
+ it->second.set_discount(d);
+ it->second.set_strength(strength);
+ }
+ }
int DecrementRule(const TRule& rule, MT19937* rng) {
RuleModelHash::iterator it = r.find(rule.f_);
@@ -46,7 +91,7 @@ struct MConditionalTranslationModel {
int IncrementRule(const TRule& rule, MT19937* rng) {
RuleModelHash::iterator it = r.find(rule.f_);
if (it == r.end()) {
- it = r.insert(make_pair(rule.f_, MFCR<1,TRule>(1.0, 1.0, 1.0, 1.0, 1e-9, 4.0))).first;
+ it = r.insert(make_pair(rule.f_, MFCR<1,TRule>(d, strength))).first;
}
p0s[0] = rp0(rule);
TableCount delta = it->second.increment(rule, p0s.begin(), lambdas.begin(), rng);
@@ -66,15 +111,7 @@ struct MConditionalTranslationModel {
}
prob_t Likelihood() const {
- prob_t p = prob_t::One();
-#if 0
- for (RuleModelHash::const_iterator it = r.begin(); it != r.end(); ++it) {
- prob_t q; q.logeq(it->second.log_crp_prob());
- p *= q;
- for (CCRP_NoTable<TRule>::const_iterator i2 = it->second.begin(); i2 != it->second.end(); ++i2)
- p *= rp0(i2->first);
- }
-#endif
+ prob_t p; p.logeq(log_likelihood(d, strength));
return p;
}
@@ -83,6 +120,7 @@ struct MConditionalTranslationModel {
MFCR<1, TRule>,
boost::hash<std::vector<WordID> > > RuleModelHash;
RuleModelHash r;
+ double d, strength;
std::vector<prob_t> lambdas;
mutable std::vector<prob_t> p0s;
};
diff --git a/gi/pf/pyp_lm.cc b/gi/pf/pyp_lm.cc
index 7ebada13..104f356b 100644
--- a/gi/pf/pyp_lm.cc
+++ b/gi/pf/pyp_lm.cc
@@ -18,7 +18,7 @@
// I use templates to handle the recursive formalation of the prior, so
// the order of the model has to be specified here, at compile time:
-#define kORDER 3
+#define kORDER 4
using namespace std;
using namespace tr1;
@@ -114,7 +114,7 @@ template <unsigned N> struct PYPLM {
if (aa <= -dd) return -std::numeric_limits<double>::infinity();
//double llh = Md::log_beta_density(dd, 10, 3) + Md::log_gamma_density(aa, 1, 1);
double llh = Md::log_beta_density(dd, discount_a, discount_b) +
- Md::log_gamma_density(aa, strength_s, strength_r);
+ Md::log_gamma_density(aa + dd, strength_s, strength_r);
typename unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > >::const_iterator it;
for (it = p.begin(); it != p.end(); ++it)
llh += it->second.log_crp_prob(dd, aa);
@@ -141,12 +141,14 @@ template <unsigned N> struct PYPLM {
DiscountResampler dr(*this);
AlphaResampler ar(*this);
for (int iter = 0; iter < nloop; ++iter) {
- strength = slice_sampler1d(ar, strength, *rng, 0.0,
+ strength = slice_sampler1d(ar, strength, *rng, -d + std::numeric_limits<double>::min(),
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
- d = slice_sampler1d(dr, d, *rng, std::numeric_limits<double>::min(),
+ double min_discount = std::numeric_limits<double>::min();
+ if (strength < 0.0) min_discount -= strength;
+ d = slice_sampler1d(dr, d, *rng, min_discount,
1.0, 0.0, niterations, 100*niterations);
}
- strength = slice_sampler1d(ar, strength, *rng, 0.0,
+ strength = slice_sampler1d(ar, strength, *rng, -d + std::numeric_limits<double>::min(),
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
typename unordered_map<vector<WordID>, CCRP<WordID>, boost::hash<vector<WordID> > >::iterator it;
cerr << "PYPLM<" << N << ">(d=" << d << ",a=" << strength << ") = " << log_likelihood(d, strength) << endl;
diff --git a/utils/ccrp.h b/utils/ccrp.h
index e24130ac..439d7e1e 100644
--- a/utils/ccrp.h
+++ b/utils/ccrp.h
@@ -225,12 +225,12 @@ class CCRP {
StrengthResampler sr(*this);
for (int iter = 0; iter < nloop; ++iter) {
if (has_strength_prior()) {
- strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
+ strength_ = slice_sampler1d(sr, strength_, *rng, -discount_ + std::numeric_limits<double>::min(),
std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
}
if (has_discount_prior()) {
double min_discount = std::numeric_limits<double>::min();
- if (strength_ < 0.0) min_discount = -strength_;
+ if (strength_ < 0.0) min_discount -= strength_;
discount_ = slice_sampler1d(dr, discount_, *rng, min_discount,
1.0, 0.0, niterations, 100*niterations);
}
diff --git a/utils/mfcr.h b/utils/mfcr.h
index 6cc0ebf1..886f01ef 100644
--- a/utils/mfcr.h
+++ b/utils/mfcr.h
@@ -48,7 +48,7 @@ class MFCR {
discount_prior_strength_(std::numeric_limits<double>::quiet_NaN()),
discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
strength_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
- strength_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+ strength_prior_rate_(std::numeric_limits<double>::quiet_NaN()) { check_hyperparameters(); }
MFCR(double discount_strength, double discount_beta, double strength_shape, double strength_rate, double d = 0.9, double strength = 10.0) :
num_tables_(),
@@ -58,10 +58,23 @@ class MFCR {
discount_prior_strength_(discount_strength),
discount_prior_beta_(discount_beta),
strength_prior_shape_(strength_shape),
- strength_prior_rate_(strength_rate) {}
+ strength_prior_rate_(strength_rate) { check_hyperparameters(); }
+
+ void check_hyperparameters() {
+ if (discount_ < 0.0 || discount_ >= 1.0) {
+ std::cerr << "Bad discount: " << discount_ << std::endl;
+ abort();
+ }
+ if (strength_ <= -discount_) {
+ std::cerr << "Bad strength: " << strength_ << " (discount=" << discount_ << ")" << std::endl;
+ abort();
+ }
+ }
double discount() const { return discount_; }
double strength() const { return strength_; }
+ void set_discount(double d) { discount_ = d; check_hyperparameters(); }
+ void set_strength(double a) { strength_ = a; check_hyperparameters(); }
bool has_discount_prior() const {
return !std::isnan(discount_prior_strength_);
@@ -275,7 +288,7 @@ class MFCR {
}
if (has_discount_prior()) {
double min_discount = std::numeric_limits<double>::min();
- if (strength_ < 0.0) min_discount = -strength_;
+ if (strength_ < 0.0) min_discount -= strength_;
discount_ = slice_sampler1d(dr, discount_, *rng, min_discount,
1.0, 0.0, niterations, 100*niterations);
}
--
cgit v1.2.3