1 #ifndef STAN_MCMC_CHAINS_HPP 2 #define STAN_MCMC_CHAINS_HPP 5 #include <stan/math/prim/mat.hpp> 6 #include <boost/accumulators/accumulators.hpp> 7 #include <boost/accumulators/statistics/stats.hpp> 8 #include <boost/accumulators/statistics/mean.hpp> 9 #include <boost/accumulators/statistics/tail_quantile.hpp> 10 #include <boost/accumulators/statistics/p_square_quantile.hpp> 11 #include <boost/accumulators/statistics/variance.hpp> 12 #include <boost/accumulators/statistics/covariance.hpp> 13 #include <boost/accumulators/statistics/variates/covariate.hpp> 14 #include <boost/random/uniform_int_distribution.hpp> 15 #include <boost/random/additive_combine.hpp> 45 template <
class RNG = boost::random::ecuyer1988>
48 Eigen::Matrix<std::string, Dynamic, 1> param_names_;
49 Eigen::Matrix<Eigen::MatrixXd, Dynamic, 1> samples_;
50 Eigen::VectorXi warmup_;
52 static double mean(
const Eigen::VectorXd& x) {
53 return (x.array() / x.size()).sum();
56 static double variance(
const Eigen::VectorXd& x) {
58 return ((x.array() - m) / std::sqrt((x.size() - 1.0))).square().sum();
61 static double sd(
const Eigen::VectorXd& x) {
62 return std::sqrt(variance(x));
66 static double covariance(
const Eigen::VectorXd& x,
67 const Eigen::VectorXd& y,
68 std::ostream* err = 0) {
69 if (x.rows() != y.rows() && err)
70 *err <<
"warning: covariance of different length chains";
71 using boost::accumulators::accumulator_set;
72 using boost::accumulators::stats;
73 using boost::accumulators::tag::variance;
74 using boost::accumulators::tag::covariance;
75 using boost::accumulators::tag::covariate1;
77 accumulator_set<double, stats<covariance<double, covariate1> > > acc;
79 int M = std::min(x.size(), y.size());
80 for (
int i = 0; i < M; i++)
81 acc(x(i), boost::accumulators::covariate1 = y(i));
83 return boost::accumulators::covariance(acc) * M / (M-1);
86 static double correlation(
const Eigen::VectorXd& x,
87 const Eigen::VectorXd& y,
88 std::ostream* err = 0) {
89 if (x.rows() != y.rows() && err)
90 *err <<
"warning: covariance of different length chains";
91 using boost::accumulators::accumulator_set;
92 using boost::accumulators::stats;
93 using boost::accumulators::tag::variance;
94 using boost::accumulators::tag::covariance;
95 using boost::accumulators::tag::covariate1;
97 accumulator_set<double, stats<variance,
98 covariance<double, covariate1> > > acc_xy;
99 accumulator_set<double, stats<variance> > acc_y;
101 int M = std::min(x.size(), y.size());
102 for (
int i = 0; i < M; i++) {
103 acc_xy(x(i), boost::accumulators::covariate1 = y(i));
107 double cov = boost::accumulators::covariance(acc_xy);
108 if (cov > -1e-8 && cov < 1e-8)
110 return cov / std::sqrt(boost::accumulators::variance(acc_xy)
111 * boost::accumulators::variance(acc_y));
114 static double quantile(
const Eigen::VectorXd& x,
const double prob) {
115 using boost::accumulators::accumulator_set;
116 using boost::accumulators::left;
117 using boost::accumulators::quantile;
118 using boost::accumulators::quantile_probability;
119 using boost::accumulators::right;
120 using boost::accumulators::stats;
121 using boost::accumulators::tag::tail;
122 using boost::accumulators::tag::tail_quantile;
125 size_t cache_size = M;
128 accumulator_set<double, stats<tail_quantile<left> > >
129 acc(tail<left>::cache_size = cache_size);
130 for (
int i = 0; i < M; i++)
132 return quantile(acc, quantile_probability = prob);
134 accumulator_set<double, stats<tail_quantile<right> > >
135 acc(tail<right>::cache_size = cache_size);
136 for (
int i = 0; i < M; i++)
138 return quantile(acc, quantile_probability = prob);
141 static Eigen::VectorXd
142 quantiles(
const Eigen::VectorXd& x,
const Eigen::VectorXd& probs) {
143 using boost::accumulators::accumulator_set;
144 using boost::accumulators::left;
145 using boost::accumulators::quantile_probability;
146 using boost::accumulators::quantile;
147 using boost::accumulators::right;
148 using boost::accumulators::stats;
149 using boost::accumulators::tag::tail;
150 using boost::accumulators::tag::tail_quantile;
154 size_t cache_size = M;
156 accumulator_set<double, stats<tail_quantile<left> > >
157 acc_left(tail<left>::cache_size = cache_size);
158 accumulator_set<double, stats<tail_quantile<right> > >
159 acc_right(tail<right>::cache_size = cache_size);
161 for (
int i = 0; i < M; i++) {
166 Eigen::VectorXd q(probs.size());
167 for (
int i = 0; i < probs.size(); i++) {
169 q(i) = quantile(acc_left,
170 quantile_probability = probs(i));
172 q(i) = quantile(acc_right,
173 quantile_probability = probs(i));
178 static Eigen::VectorXd autocorrelation(
const Eigen::VectorXd& x) {
180 using stan::math::index_type;
181 typedef typename index_type<vector<double> >::type idx_t;
183 std::vector<double> ac;
184 std::vector<double>
sample(x.size());
185 for (
int i = 0; i < x.size(); i++)
187 stan::math::autocorrelation(
sample, ac);
189 Eigen::VectorXd ac2(ac.size());
190 for (idx_t i = 0; i < ac.size(); i++)
195 static Eigen::VectorXd autocovariance(
const Eigen::VectorXd& x) {
197 using stan::math::index_type;
198 typedef typename index_type<vector<double> >::type idx_t;
200 std::vector<double> ac;
201 std::vector<double>
sample(x.size());
202 for (
int i = 0; i < x.size(); i++)
204 stan::math::autocovariance(
sample, ac);
206 Eigen::VectorXd ac2(ac.size());
207 for (idx_t i = 0; i < ac.size(); i++)
227 double effective_sample_size(
const Eigen::Matrix<Eigen::VectorXd,
232 int n_samples =
samples(0).size();
233 for (
int chain = 1; chain <
chains; chain++) {
234 n_samples = std::min(n_samples,
235 static_cast<int>(
samples(chain).size()));
238 Eigen::Matrix<Eigen::VectorXd, Dynamic, 1> acov(chains);
239 for (
int chain = 0; chain <
chains; chain++) {
240 acov(chain) = autocovariance(
samples(chain));
243 Eigen::VectorXd chain_mean(chains);
244 Eigen::VectorXd chain_var(chains);
245 for (
int chain = 0; chain <
chains; chain++) {
247 chain_mean(chain) = mean(
samples(chain));
248 chain_var(chain) = acov(chain)(0)*n_kept_samples/(n_kept_samples-1);
251 double mean_var = mean(chain_var);
252 double var_plus = mean_var*(n_samples-1)/n_samples;
254 var_plus += variance(chain_mean);
255 Eigen::VectorXd rho_hat_t(n_samples);
259 for (
int t = 1; (t < n_samples && rho_hat >= 0); t++) {
260 Eigen::VectorXd acov_t(chains);
261 for (
int chain = 0; chain <
chains; chain++) {
262 acov_t(chain) = acov(chain)(t);
264 rho_hat = 1 - (mean_var - mean(acov_t)) / var_plus;
266 rho_hat_t(t) = rho_hat;
269 double ess = chains * n_samples;
271 ess /= 1 + 2 * rho_hat_t.sum();
290 split_potential_scale_reduction(
291 const Eigen::Matrix<Eigen::VectorXd,
294 int n_samples =
samples(0).size();
295 for (
int chain = 1; chain <
chains; chain++) {
296 n_samples = std::min(n_samples,
297 static_cast<int>(
samples(chain).size()));
299 if (n_samples % 2 == 1)
301 int n = n_samples / 2;
303 Eigen::VectorXd split_chain_mean(2*chains);
304 Eigen::VectorXd split_chain_var(2*chains);
306 for (
int chain = 0; chain <
chains; chain++) {
307 split_chain_mean(2*chain) = mean(
samples(chain).topRows(n));
308 split_chain_mean(2*chain+1) = mean(
samples(chain).bottomRows(n));
310 split_chain_var(2*chain) = variance(
samples(chain).topRows(n));
311 split_chain_var(2*chain+1) = variance(
samples(chain).bottomRows(n));
314 double var_between = n * variance(split_chain_mean);
315 double var_within = mean(split_chain_var);
318 return sqrt((var_between/var_within + n-1)/n);
323 : param_names_(param_names) { }
326 : param_names_(param_names.size()) {
327 for (
size_t i = 0; i < param_names.size(); i++)
328 param_names_(i) = param_names[i];
332 : param_names_(stan_csv.header) {
333 if (stan_csv.
samples.rows() > 0)
338 return samples_.size();
342 return param_names_.size();
345 const Eigen::Matrix<std::string, Dynamic, 1>&
param_names()
const {
350 return param_names_(j);
353 int index(
const std::string& name)
const {
355 for (
int i = 0; i < param_names_.size(); i++)
356 if (param_names_(i) == name)
366 warmup_.setConstant(warmup);
374 return warmup_(chain);
378 return samples_(chain).rows();
383 for (
int chain = 0; chain <
num_chains(); chain++)
394 for (
int chain = 0; chain <
num_chains(); chain++)
400 const Eigen::MatrixXd&
sample) {
402 throw std::invalid_argument(
"add(chain, sample): number of columns" 403 " in sample does not match chains");
409 Eigen::Matrix<Eigen::MatrixXd, Dynamic, 1>
412 for (
int i = 0; i < n; i++) {
413 samples_copy(i) = samples_(i);
414 warmup_copy(i) = warmup_(i);
417 samples_.resize(chain+1);
418 warmup_.resize(chain+1);
419 for (
int i = 0; i < n; i++) {
420 samples_(i) = samples_copy(i);
421 warmup_(i) = warmup_copy(i);
423 for (
int i = n; i < chain+1; i++) {
424 samples_(i) = Eigen::MatrixXd(0,
num_params());
428 int row = samples_(chain).rows();
429 Eigen::MatrixXd new_samples(row+sample.rows(),
num_params());
430 new_samples << samples_(chain),
sample;
431 samples_(chain) = new_samples;
435 if (sample.rows() == 0)
438 throw std::invalid_argument(
"add(sample): number of columns in" 439 " sample does not match chains");
450 void add(
const std::vector<std::vector<double> >&
sample) {
451 int n_row =
sample.size();
454 int n_col =
sample[0].size();
455 Eigen::MatrixXd sample_copy(n_row, n_col);
456 for (
int i = 0; i < n_row; i++) {
458 = Eigen::VectorXd::Map(&
sample[i][0],
sample[0].size());
465 throw std::invalid_argument(
"add(stan_csv): number of columns in" 466 " sample does not match chains");
467 if (!param_names_.cwiseEqual(stan_csv.
header).all()) {
468 throw std::invalid_argument(
"add(stan_csv): header does not match" 483 for (
int chain = 0; chain <
num_chains(); chain++) {
485 s.middleRows(start, n) = samples_(chain).col(index).bottomRows(n);
491 Eigen::VectorXd
samples(
const int chain,
const std::string& name)
const {
495 Eigen::VectorXd
samples(
const std::string& name)
const {
500 return mean(
samples(chain, index));
507 double mean(
const int chain,
const std::string& name)
const {
508 return mean(chain,
index(name));
511 double mean(
const std::string& name)
const {
512 return mean(
index(name));
515 double sd(
const int chain,
const int index)
const {
516 return sd(
samples(chain, index));
523 double sd(
const int chain,
const std::string& name)
const {
524 return sd(chain,
index(name));
527 double sd(
const std::string& name)
const {
528 return sd(
index(name));
532 return variance(
samples(chain, index));
536 return variance(
samples(index));
539 double variance(
const int chain,
const std::string& name)
const {
540 return variance(chain,
index(name));
544 return variance(
index(name));
548 covariance(
const int chain,
const int index1,
const int index2)
const {
549 return covariance(
samples(chain, index1),
samples(chain, index2));
552 double covariance(
const int index1,
const int index2)
const {
557 const std::string& name2)
const {
558 return covariance(chain,
index(name1),
index(name2));
562 covariance(
const std::string& name1,
const std::string& name2)
const {
563 return covariance(
index(name1),
index(name2));
567 correlation(
const int chain,
const int index1,
const int index2)
const {
568 return correlation(
samples(chain, index1),
samples(chain, index2));
576 const std::string& name2)
const {
577 return correlation(chain,
index(name1),
index(name2));
581 correlation(
const std::string& name1,
const std::string& name2)
const {
582 return correlation(
index(name1),
index(name2));
587 return quantile(
samples(chain, index), prob);
591 return quantile(
samples(index), prob);
594 double quantile(
int chain,
const std::string& name,
double prob)
const {
595 return quantile(chain,
index(name), prob);
598 double quantile(
const std::string& name,
const double prob)
const {
599 return quantile(
index(name), prob);
604 return quantiles(
samples(chain, index), probs);
608 return quantiles(
samples(index), probs);
613 const Eigen::VectorXd& probs)
const {
614 return quantiles(chain,
index(name), probs);
618 quantiles(
const std::string& name,
const Eigen::VectorXd& probs)
const {
619 return quantiles(
index(name), probs);
624 double low_prob = (1-prob)/2;
625 double high_prob = 1-low_prob;
627 Eigen::Vector2d interval;
629 << quantile(chain, index, low_prob),
630 quantile(chain, index, high_prob);
635 double low_prob = (1-prob)/2;
636 double high_prob = 1-low_prob;
638 Eigen::Vector2d interval;
639 interval << quantile(index, low_prob), quantile(index, high_prob);
655 return autocorrelation(
samples(chain, index));
659 const std::string& name)
const {
660 return autocorrelation(chain,
index(name));
664 return autocovariance(
samples(chain, index));
668 return autocovariance(chain,
index(name));
673 Eigen::Matrix<Eigen::VectorXd, Dynamic, 1>
675 for (
int chain = 0; chain <
num_chains(); chain++) {
678 return effective_sample_size(samples);
682 return effective_sample_size(
index(name));
686 Eigen::Matrix<Eigen::VectorXd, Dynamic, 1>
688 for (
int chain = 0; chain <
num_chains(); chain++) {
691 return split_potential_scale_reduction(samples);
695 return split_potential_scale_reduction(
index(name));
double covariance(const int chain, const int index1, const int index2) const
double variance(const std::string &name) const
int num_kept_samples() const
double effective_sample_size(const std::string &name) const
int index(const std::string &name) const
double correlation(const int chain, const int index1, const int index2) const
double correlation(const int chain, const std::string &name1, const std::string &name2) const
int num_kept_samples(const int chain) const
void set_warmup(const int warmup)
const Eigen::VectorXi & warmup() const
double covariance(const int chain, const std::string &name1, const std::string &name2) const
Eigen::VectorXd samples(const std::string &name) const
Eigen::VectorXd autocorrelation(const int chain, const int index) const
double quantile(const int chain, const int index, const double prob) const
void sample(stan::mcmc::base_mcmc *sampler, int num_warmup, int num_samples, int num_thin, int refresh, bool save, stan::services::sample::mcmc_writer< Model, SampleRecorder, DiagnosticRecorder, MessageRecorder > &mcmc_writer, stan::mcmc::sample &init_s, Model &model, RNG &base_rng, const std::string &prefix, const std::string &suffix, std::ostream &o, StartTransitionCallback &callback, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Probability, optimization and sampling library.
Eigen::VectorXd samples(const int chain, const std::string &name) const
int num_samples(const int chain) const
Eigen::VectorXd quantiles(int chain, int index, const Eigen::VectorXd &probs) const
chains(const stan::io::stan_csv &stan_csv)
Eigen::Vector2d central_interval(const std::string &name, double prob) const
void add(const stan::io::stan_csv &stan_csv)
Eigen::VectorXd quantiles(int chain, const std::string &name, const Eigen::VectorXd &probs) const
double mean(const int index) const
double sd(const int index) const
double correlation(const std::string &name1, const std::string &name2) const
Eigen::Vector2d central_interval(int chain, const std::string &name, double prob) const
Eigen::VectorXd autocorrelation(int chain, const std::string &name) const
double covariance(const std::string &name1, const std::string &name2) const
double split_potential_scale_reduction(const int index) const
Eigen::VectorXd samples(const int chain, const int index) const
Eigen::VectorXd samples(const int index) const
chains(const Eigen::Matrix< std::string, Dynamic, 1 > ¶m_names)
double sd(const std::string &name) const
const Eigen::Matrix< std::string, Dynamic, 1 > & param_names() const
double correlation(const int index1, const int index2) const
double variance(const int index) const
void add(const std::vector< std::vector< double > > &sample)
Convert a vector of vector<double> to Eigen::MatrixXd.
const std::string & param_name(int j) const
Eigen::VectorXd quantiles(const std::string &name, const Eigen::VectorXd &probs) const
double sd(const int chain, const int index) const
double quantile(int chain, const std::string &name, double prob) const
chains(const std::vector< std::string > ¶m_names)
int warmup(const int chain) const
Eigen::VectorXd autocovariance(int chain, const std::string &name) const
double variance(const int chain, const int index) const
double effective_sample_size(const int index) const
void add(const Eigen::MatrixXd &sample)
An mcmc::chains object stores parameter names and dimensionalities along with samples from multiple c...
double covariance(const int index1, const int index2) const
double quantile(const int index, const double prob) const
double split_potential_scale_reduction(const std::string &name) const
void add(const int chain, const Eigen::MatrixXd &sample)
double mean(const std::string &name) const
Eigen::VectorXd autocovariance(const int chain, const int index) const
Eigen::Vector2d central_interval(int chain, int index, double prob) const
double mean(const int chain, const std::string &name) const
Eigen::Vector2d central_interval(int index, double prob) const
double variance(const int chain, const std::string &name) const
Eigen::VectorXd quantiles(int index, const Eigen::VectorXd &probs) const
double sd(const int chain, const std::string &name) const
stan_csv_metadata metadata
double mean(const int chain, const int index) const
Eigen::Matrix< std::string, Eigen::Dynamic, 1 > header
double quantile(const std::string &name, const double prob) const
void set_warmup(const int chain, const int warmup)