1 #ifndef STAN_MCMC_CHAINS_HPP
2 #define STAN_MCMC_CHAINS_HPP
5 #include <stan/math/prim/mat/fun/variance.hpp>
6 #include <stan/math/prim/arr/meta/index_type.hpp>
7 #include <stan/math/prim/mat/meta/index_type.hpp>
8 #include <stan/math/prim/scal/meta/index_type.hpp>
9 #include <stan/math/prim/mat/fun/autocorrelation.hpp>
10 #include <stan/math/prim/mat/fun/autocovariance.hpp>
11 #include <boost/accumulators/accumulators.hpp>
12 #include <boost/accumulators/statistics/stats.hpp>
13 #include <boost/accumulators/statistics/mean.hpp>
14 #include <boost/accumulators/statistics/tail_quantile.hpp>
15 #include <boost/accumulators/statistics/p_square_quantile.hpp>
16 #include <boost/accumulators/statistics/variance.hpp>
17 #include <boost/accumulators/statistics/covariance.hpp>
18 #include <boost/accumulators/statistics/variates/covariate.hpp>
19 #include <boost/random/uniform_int_distribution.hpp>
20 #include <boost/random/additive_combine.hpp>
50 template <
class RNG = boost::random::ecuyer1988>
53 Eigen::Matrix<std::string, Dynamic, 1> param_names_;
54 Eigen::Matrix<Eigen::MatrixXd, Dynamic, 1> samples_;
55 Eigen::VectorXi warmup_;
57 static double mean(
const Eigen::VectorXd& x) {
58 return (x.array() / x.size()).sum();
61 static double variance(
const Eigen::VectorXd& x) {
63 return ((x.array() - m) / std::sqrt((x.size() - 1.0))).square().sum();
66 static double sd(
const Eigen::VectorXd& x) {
67 return std::sqrt(variance(x));
71 static double covariance(
const Eigen::VectorXd& x,
72 const Eigen::VectorXd& y,
73 std::ostream* err = 0) {
74 if (x.rows() != y.rows() && err)
75 *err <<
"warning: covariance of different length chains";
76 using boost::accumulators::accumulator_set;
77 using boost::accumulators::stats;
78 using boost::accumulators::tag::variance;
79 using boost::accumulators::tag::covariance;
80 using boost::accumulators::tag::covariate1;
82 accumulator_set<double, stats<covariance<double, covariate1> > > acc;
84 int M = std::min(x.size(), y.size());
85 for (
int i = 0; i < M; i++)
86 acc(x(i), boost::accumulators::covariate1 = y(i));
88 return boost::accumulators::covariance(acc) * M / (M-1);
91 static double correlation(
const Eigen::VectorXd& x,
92 const Eigen::VectorXd& y,
93 std::ostream* err = 0) {
94 if (x.rows() != y.rows() && err)
95 *err <<
"warning: covariance of different length chains";
96 using boost::accumulators::accumulator_set;
97 using boost::accumulators::stats;
98 using boost::accumulators::tag::variance;
99 using boost::accumulators::tag::covariance;
100 using boost::accumulators::tag::covariate1;
102 accumulator_set<double, stats<variance,
103 covariance<double, covariate1> > > acc_xy;
104 accumulator_set<double, stats<variance> > acc_y;
106 int M = std::min(x.size(), y.size());
107 for (
int i = 0; i < M; i++) {
108 acc_xy(x(i), boost::accumulators::covariate1 = y(i));
112 double cov = boost::accumulators::covariance(acc_xy);
113 if (cov > -1e-8 && cov < 1e-8)
115 return cov / std::sqrt(boost::accumulators::variance(acc_xy)
116 * boost::accumulators::variance(acc_y));
119 static double quantile(
const Eigen::VectorXd& x,
const double prob) {
120 using boost::accumulators::accumulator_set;
121 using boost::accumulators::left;
122 using boost::accumulators::quantile;
123 using boost::accumulators::quantile_probability;
124 using boost::accumulators::right;
125 using boost::accumulators::stats;
126 using boost::accumulators::tag::tail;
127 using boost::accumulators::tag::tail_quantile;
130 size_t cache_size = M;
133 accumulator_set<double, stats<tail_quantile<left> > >
134 acc(tail<left>::cache_size = cache_size);
135 for (
int i = 0; i < M; i++)
137 return quantile(acc, quantile_probability = prob);
139 accumulator_set<double, stats<tail_quantile<right> > >
140 acc(tail<right>::cache_size = cache_size);
141 for (
int i = 0; i < M; i++)
143 return quantile(acc, quantile_probability = prob);
146 static Eigen::VectorXd
147 quantiles(
const Eigen::VectorXd& x,
const Eigen::VectorXd& probs) {
148 using boost::accumulators::accumulator_set;
149 using boost::accumulators::left;
150 using boost::accumulators::quantile_probability;
151 using boost::accumulators::quantile;
152 using boost::accumulators::right;
153 using boost::accumulators::stats;
154 using boost::accumulators::tag::tail;
155 using boost::accumulators::tag::tail_quantile;
159 size_t cache_size = M;
161 accumulator_set<double, stats<tail_quantile<left> > >
162 acc_left(tail<left>::cache_size = cache_size);
163 accumulator_set<double, stats<tail_quantile<right> > >
164 acc_right(tail<right>::cache_size = cache_size);
166 for (
int i = 0; i < M; i++) {
171 Eigen::VectorXd q(probs.size());
172 for (
int i = 0; i < probs.size(); i++) {
174 q(i) = quantile(acc_left,
175 quantile_probability = probs(i));
177 q(i) = quantile(acc_right,
178 quantile_probability = probs(i));
183 static Eigen::VectorXd autocorrelation(
const Eigen::VectorXd& x) {
185 using stan::math::index_type;
186 typedef typename index_type<vector<double> >::type idx_t;
188 std::vector<double> ac;
189 std::vector<double>
sample(x.size());
190 for (
int i = 0; i < x.size(); i++)
192 stan::math::autocorrelation(
sample, ac);
194 Eigen::VectorXd ac2(ac.size());
195 for (idx_t i = 0; i < ac.size(); i++)
200 static Eigen::VectorXd autocovariance(
const Eigen::VectorXd& x) {
202 using stan::math::index_type;
203 typedef typename index_type<vector<double> >::type idx_t;
205 std::vector<double> ac;
206 std::vector<double>
sample(x.size());
207 for (
int i = 0; i < x.size(); i++)
209 stan::math::autocovariance(
sample, ac);
211 Eigen::VectorXd ac2(ac.size());
212 for (idx_t i = 0; i < ac.size(); i++)
232 double effective_sample_size(
const Eigen::Matrix<Eigen::VectorXd,
237 int n_samples =
samples(0).size();
238 for (
int chain = 1; chain <
chains; chain++) {
239 n_samples = std::min(n_samples,
240 static_cast<int>(
samples(chain).size()));
243 Eigen::Matrix<Eigen::VectorXd, Dynamic, 1> acov(chains);
244 for (
int chain = 0; chain <
chains; chain++) {
245 acov(chain) = autocovariance(
samples(chain));
248 Eigen::VectorXd chain_mean(chains);
249 Eigen::VectorXd chain_var(chains);
250 for (
int chain = 0; chain <
chains; chain++) {
252 chain_mean(chain) = mean(
samples(chain));
253 chain_var(chain) = acov(chain)(0)*n_kept_samples/(n_kept_samples-1);
256 double mean_var = mean(chain_var);
257 double var_plus = mean_var*(n_samples-1)/n_samples;
259 var_plus += variance(chain_mean);
260 Eigen::VectorXd rho_hat_t(n_samples);
264 for (
int t = 1; (t < n_samples && rho_hat >= 0); t++) {
265 Eigen::VectorXd acov_t(chains);
266 for (
int chain = 0; chain <
chains; chain++) {
267 acov_t(chain) = acov(chain)(t);
269 rho_hat = 1 - (mean_var - mean(acov_t)) / var_plus;
271 rho_hat_t(t) = rho_hat;
274 double ess = chains * n_samples;
276 ess /= 1 + 2 * rho_hat_t.sum();
295 split_potential_scale_reduction(
296 const Eigen::Matrix<Eigen::VectorXd,
299 int n_samples =
samples(0).size();
300 for (
int chain = 1; chain <
chains; chain++) {
301 n_samples = std::min(n_samples,
302 static_cast<int>(
samples(chain).size()));
304 if (n_samples % 2 == 1)
306 int n = n_samples / 2;
308 Eigen::VectorXd split_chain_mean(2*chains);
309 Eigen::VectorXd split_chain_var(2*chains);
311 for (
int chain = 0; chain <
chains; chain++) {
312 split_chain_mean(2*chain) = mean(
samples(chain).topRows(n));
313 split_chain_mean(2*chain+1) = mean(
samples(chain).bottomRows(n));
315 split_chain_var(2*chain) = variance(
samples(chain).topRows(n));
316 split_chain_var(2*chain+1) = variance(
samples(chain).bottomRows(n));
319 double var_between = n * variance(split_chain_mean);
320 double var_within = mean(split_chain_var);
323 return sqrt((var_between/var_within + n-1)/n);
328 : param_names_(param_names) { }
331 : param_names_(param_names.size()) {
332 for (
size_t i = 0; i < param_names.size(); i++)
333 param_names_(i) = param_names[i];
337 : param_names_(stan_csv.header) {
338 if (stan_csv.
samples.rows() > 0)
343 return samples_.size();
347 return param_names_.size();
350 const Eigen::Matrix<std::string, Dynamic, 1>&
param_names()
const {
355 return param_names_(j);
358 int index(
const std::string& name)
const {
360 for (
int i = 0; i < param_names_.size(); i++)
361 if (param_names_(i) == name)
371 warmup_.setConstant(warmup);
379 return warmup_(chain);
383 return samples_(chain).rows();
388 for (
int chain = 0; chain <
num_chains(); chain++)
399 for (
int chain = 0; chain <
num_chains(); chain++)
405 const Eigen::MatrixXd&
sample) {
407 throw std::invalid_argument(
"add(chain, sample): number of columns"
408 " in sample does not match chains");
414 Eigen::Matrix<Eigen::MatrixXd, Dynamic, 1>
417 for (
int i = 0; i < n; i++) {
418 samples_copy(i) = samples_(i);
419 warmup_copy(i) = warmup_(i);
422 samples_.resize(chain+1);
423 warmup_.resize(chain+1);
424 for (
int i = 0; i < n; i++) {
425 samples_(i) = samples_copy(i);
426 warmup_(i) = warmup_copy(i);
428 for (
int i = n; i < chain+1; i++) {
429 samples_(i) = Eigen::MatrixXd(0,
num_params());
433 int row = samples_(chain).rows();
434 Eigen::MatrixXd new_samples(row+sample.rows(),
num_params());
435 new_samples << samples_(chain),
sample;
436 samples_(chain) = new_samples;
440 if (sample.rows() == 0)
443 throw std::invalid_argument(
"add(sample): number of columns in"
444 " sample does not match chains");
455 void add(
const std::vector<std::vector<double> >&
sample) {
456 int n_row =
sample.size();
459 int n_col =
sample[0].size();
460 Eigen::MatrixXd sample_copy(n_row, n_col);
461 for (
int i = 0; i < n_row; i++) {
463 = Eigen::VectorXd::Map(&
sample[i][0],
sample[0].size());
470 throw std::invalid_argument(
"add(stan_csv): number of columns in"
471 " sample does not match chains");
472 if (!param_names_.cwiseEqual(stan_csv.
header).all()) {
473 throw std::invalid_argument(
"add(stan_csv): header does not match"
488 for (
int chain = 0; chain <
num_chains(); chain++) {
490 s.middleRows(start, n) = samples_(chain).col(index).bottomRows(n);
496 Eigen::VectorXd
samples(
const int chain,
const std::string& name)
const {
500 Eigen::VectorXd
samples(
const std::string& name)
const {
505 return mean(
samples(chain, index));
512 double mean(
const int chain,
const std::string& name)
const {
513 return mean(chain,
index(name));
516 double mean(
const std::string& name)
const {
517 return mean(
index(name));
520 double sd(
const int chain,
const int index)
const {
521 return sd(
samples(chain, index));
528 double sd(
const int chain,
const std::string& name)
const {
529 return sd(chain,
index(name));
532 double sd(
const std::string& name)
const {
533 return sd(
index(name));
537 return variance(
samples(chain, index));
541 return variance(
samples(index));
544 double variance(
const int chain,
const std::string& name)
const {
545 return variance(chain,
index(name));
549 return variance(
index(name));
553 covariance(
const int chain,
const int index1,
const int index2)
const {
554 return covariance(
samples(chain, index1),
samples(chain, index2));
557 double covariance(
const int index1,
const int index2)
const {
562 const std::string& name2)
const {
563 return covariance(chain,
index(name1),
index(name2));
567 covariance(
const std::string& name1,
const std::string& name2)
const {
568 return covariance(
index(name1),
index(name2));
572 correlation(
const int chain,
const int index1,
const int index2)
const {
573 return correlation(
samples(chain, index1),
samples(chain, index2));
581 const std::string& name2)
const {
582 return correlation(chain,
index(name1),
index(name2));
586 correlation(
const std::string& name1,
const std::string& name2)
const {
587 return correlation(
index(name1),
index(name2));
592 return quantile(
samples(chain, index), prob);
596 return quantile(
samples(index), prob);
599 double quantile(
int chain,
const std::string& name,
double prob)
const {
600 return quantile(chain,
index(name), prob);
603 double quantile(
const std::string& name,
const double prob)
const {
604 return quantile(
index(name), prob);
609 return quantiles(
samples(chain, index), probs);
613 return quantiles(
samples(index), probs);
618 const Eigen::VectorXd& probs)
const {
619 return quantiles(chain,
index(name), probs);
623 quantiles(
const std::string& name,
const Eigen::VectorXd& probs)
const {
624 return quantiles(
index(name), probs);
629 double low_prob = (1-prob)/2;
630 double high_prob = 1-low_prob;
632 Eigen::Vector2d interval;
634 << quantile(chain, index, low_prob),
635 quantile(chain, index, high_prob);
640 double low_prob = (1-prob)/2;
641 double high_prob = 1-low_prob;
643 Eigen::Vector2d interval;
644 interval << quantile(index, low_prob), quantile(index, high_prob);
660 return autocorrelation(
samples(chain, index));
664 const std::string& name)
const {
665 return autocorrelation(chain,
index(name));
669 return autocovariance(
samples(chain, index));
673 return autocovariance(chain,
index(name));
678 Eigen::Matrix<Eigen::VectorXd, Dynamic, 1>
680 for (
int chain = 0; chain <
num_chains(); chain++) {
683 return effective_sample_size(samples);
687 return effective_sample_size(
index(name));
691 Eigen::Matrix<Eigen::VectorXd, Dynamic, 1>
693 for (
int chain = 0; chain <
num_chains(); chain++) {
696 return split_potential_scale_reduction(samples);
700 return split_potential_scale_reduction(
index(name));
int warmup(const int chain) const
int num_kept_samples() const
double sd(const std::string &name) const
void set_warmup(const int warmup)
Eigen::VectorXd samples(const int chain, const std::string &name) const
double covariance(const std::string &name1, const std::string &name2) const
const std::string & param_name(int j) const
double correlation(const int index1, const int index2) const
double quantile(const std::string &name, const double prob) const
double mean(const int chain, const int index) const
Eigen::VectorXd quantiles(int chain, const std::string &name, const Eigen::VectorXd &probs) 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)
double quantile(const int index, const double prob) const
Probability, optimization and sampling library.
chains(const stan::io::stan_csv &stan_csv)
void add(const stan::io::stan_csv &stan_csv)
Eigen::VectorXd autocovariance(int chain, const std::string &name) const
Eigen::VectorXd samples(const int chain, const int index) const
double sd(const int chain, const int index) const
int index(const std::string &name) const
Eigen::Vector2d central_interval(int index, double prob) const
double correlation(const std::string &name1, const std::string &name2) const
double mean(const std::string &name) const
double variance(const std::string &name) const
double sd(const int index) const
double split_potential_scale_reduction(const int index) const
double variance(const int chain, const int index) const
double quantile(int chain, const std::string &name, double prob) const
double covariance(const int chain, const int index1, const int index2) const
double effective_sample_size(const int index) const
chains(const Eigen::Matrix< std::string, Dynamic, 1 > ¶m_names)
double correlation(const int chain, const int index1, const int index2) const
double mean(const int chain, const std::string &name) const
double sd(const int chain, const std::string &name) const
double covariance(const int index1, const int index2) const
void add(const std::vector< std::vector< double > > &sample)
Convert a vector of vector to Eigen::MatrixXd.
Eigen::VectorXd autocorrelation(const int chain, const int index) const
double variance(const int chain, const std::string &name) const
chains(const std::vector< std::string > ¶m_names)
Eigen::VectorXd samples(const int index) const
const Eigen::Matrix< std::string, Dynamic, 1 > & param_names() const
int num_samples(const int chain) const
double correlation(const int chain, const std::string &name1, const std::string &name2) const
void add(const Eigen::MatrixXd &sample)
double effective_sample_size(const std::string &name) const
double variance(const int index) const
double mean(const int index) const
An mcmc::chains object stores parameter names and dimensionalities along with samples from multiple c...
Eigen::Vector2d central_interval(int chain, const std::string &name, double prob) const
Eigen::VectorXd autocovariance(const int chain, const int index) const
double split_potential_scale_reduction(const std::string &name) const
void add(const int chain, const Eigen::MatrixXd &sample)
Eigen::Vector2d central_interval(const std::string &name, double prob) const
Eigen::VectorXd quantiles(int index, const Eigen::VectorXd &probs) const
Eigen::VectorXd autocorrelation(int chain, const std::string &name) const
Eigen::VectorXd quantiles(int chain, int index, const Eigen::VectorXd &probs) const
Eigen::VectorXd quantiles(const std::string &name, const Eigen::VectorXd &probs) const
stan_csv_metadata metadata
double covariance(const int chain, const std::string &name1, const std::string &name2) const
double quantile(const int chain, const int index, const double prob) const
Eigen::Matrix< std::string, Eigen::Dynamic, 1 > header
Eigen::VectorXd samples(const std::string &name) const
Eigen::Vector2d central_interval(int chain, int index, double prob) const
const Eigen::VectorXi & warmup() const
int num_kept_samples(const int chain) const
void set_warmup(const int chain, const int warmup)