3 #include <gtest/gtest.h> 4 #include <boost/random/additive_combine.hpp> 16 blocker1_stream.open(
"src/test/unit/mcmc/test_csv_files/blocker.1.csv");
17 blocker2_stream.open(
"src/test/unit/mcmc/test_csv_files/blocker.2.csv");
18 epil1_stream.open(
"src/test/unit/mcmc/test_csv_files/epil.1.csv");
19 epil2_stream.open(
"src/test/unit/mcmc/test_csv_files/epil.2.csv");
32 std::stringstream
out;
34 EXPECT_EQ(
"", out.str());
37 EXPECT_EQ(0, chains1.num_chains());
38 EXPECT_EQ(blocker1.
header.size(), chains1.num_params());
39 for (
int i = 0;
i < blocker1.
header.size();
i++)
40 EXPECT_EQ(blocker1.
header(
i), chains1.param_name(
i));
47 for (
int i = 0;
i < blocker1.
header.size();
i++)
49 EXPECT_EQ(0, chains2.
warmup(0));
55 std::stringstream
out;
57 EXPECT_EQ(
"", out.str());
61 EXPECT_EQ(0, chains.num_chains());
62 EXPECT_EQ(0, chains.num_samples());
65 EXPECT_NO_THROW(chains.add(1, theta))
66 <<
"adding a single sample to a new chain";
67 EXPECT_EQ(2, chains.num_chains());
68 EXPECT_EQ(0, chains.num_samples(0));
69 EXPECT_EQ(1, chains.num_samples(1));
70 EXPECT_EQ(1, chains.num_samples());
72 theta = blocker1.
samples.row(1);
73 EXPECT_NO_THROW(chains.add(1, theta))
74 <<
"adding a single sample to an existing chain";
75 EXPECT_EQ(2, chains.num_chains());
76 EXPECT_EQ(0, chains.num_samples(0));
77 EXPECT_EQ(2, chains.num_samples(1));
78 EXPECT_EQ(2, chains.num_samples());
80 EXPECT_NO_THROW(chains.add(3, blocker1.
samples))
81 <<
"adding multiple samples to a new chain";
82 EXPECT_EQ(4, chains.num_chains());
83 EXPECT_EQ(0, chains.num_samples(0));
84 EXPECT_EQ(2, chains.num_samples(1));
85 EXPECT_EQ(0, chains.num_samples(2));
86 EXPECT_EQ(1000, chains.num_samples(3));
87 EXPECT_EQ(1002, chains.num_samples());
89 EXPECT_NO_THROW(chains.add(3, blocker1.
samples))
90 <<
"adding multiple samples to an existing chain";
91 EXPECT_EQ(4, chains.num_chains());
92 EXPECT_EQ(0, chains.num_samples(0));
93 EXPECT_EQ(2, chains.num_samples(1));
94 EXPECT_EQ(0, chains.num_samples(2));
95 EXPECT_EQ(2000, chains.num_samples(3));
96 EXPECT_EQ(2002, chains.num_samples());
99 EXPECT_NO_THROW(chains.add(blocker1.
samples))
100 <<
"adding multiple samples, adds new chain";
101 EXPECT_EQ(5, chains.num_chains());
102 EXPECT_EQ(0, chains.num_samples(0));
103 EXPECT_EQ(2, chains.num_samples(1));
104 EXPECT_EQ(0, chains.num_samples(2));
105 EXPECT_EQ(2000, chains.num_samples(3));
106 EXPECT_EQ(1000, chains.num_samples(4));
107 EXPECT_EQ(3002, chains.num_samples());
111 EXPECT_EQ(
"", out.str());
112 theta.resize(epil1.
samples.cols());
115 <<
"adding mismatched sample to an existing chain";
117 <<
"adding mismatched sample to a new chain";
119 <<
"adding mismatched samples to an existing chain";
121 <<
"adding mismatched samples to a new chain";
123 <<
"adding mismatched sample";
125 EXPECT_EQ(5, chains.num_chains())
126 <<
"validate state is identical to before";
127 EXPECT_EQ(0, chains.num_samples(0))
128 <<
"validate state is identical to before";
129 EXPECT_EQ(2, chains.num_samples(1))
130 <<
"validate state is identical to before";
131 EXPECT_EQ(0, chains.num_samples(2))
132 <<
"validate state is identical to before";
133 EXPECT_EQ(2000, chains.num_samples(3))
134 <<
"validate state is identical to before";
135 EXPECT_EQ(1000, chains.num_samples(4))
136 <<
"validate state is identical to before";
137 EXPECT_EQ(3002, chains.num_samples())
138 <<
"validate state is identical to before";
142 std::stringstream
out;
144 EXPECT_EQ(
"", out.str());
147 std::vector<std::string> param_names(blocker1.
header.size());
148 for (
int i = 0;
i < blocker1.
header.size();
i++) {
149 param_names[
i] = blocker1.
header[
i];
156 std::vector<std::vector<double> > samples(blocker1.
samples.rows());
157 for (
int i = 0;
i < blocker1.
samples.rows();
i++) {
158 samples[
i] = std::vector<double>(blocker1.
samples.cols());
160 for (
int i = 0;
i < blocker1.
samples.rows();
i++) {
161 for (
int j = 0;
j < blocker1.
samples.cols();
j++) {
166 EXPECT_EQ(blocker1.
samples.rows(),
static_cast<int>(samples.size()));
167 EXPECT_EQ(blocker1.
samples.cols(),
static_cast<int>(samples[0].size()));
169 EXPECT_NO_THROW(chains.
add(samples))
170 <<
"adding multiple samples, adds new chain";
178 std::stringstream
out;
180 EXPECT_EQ(
"", out.str());
184 EXPECT_EQ(1, chains.num_chains());
188 std::stringstream
out;
190 EXPECT_EQ(
"", out.str());
194 EXPECT_EQ(1000, chains.num_samples());
198 std::stringstream
out;
200 EXPECT_EQ(
"", out.str());
204 EXPECT_EQ(1000, chains.num_samples());
208 std::stringstream
out;
211 EXPECT_EQ(
"", out.str());
214 chains.add(blocker2);
216 EXPECT_EQ(2000, chains.num_samples());
217 EXPECT_EQ(1000, chains.num_samples(0));
218 EXPECT_EQ(1000, chains.num_samples(1));
223 std::stringstream
out;
225 EXPECT_EQ(
"", out.str());
228 ASSERT_EQ(blocker1.
header.size(), chains.num_params());
229 ASSERT_EQ(blocker1.
header.size(), chains.param_names().size());
230 for (
int i = 0;
i < blocker1.
header.size();
i++) {
231 EXPECT_EQ(blocker1.
header(
i), chains.param_names()(
i));
235 std::stringstream
out;
237 EXPECT_EQ(
"", out.str());
240 ASSERT_EQ(blocker1.
header.size(), chains.num_params());
241 for (
int i = 0;
i < blocker1.
header.size();
i++) {
242 EXPECT_EQ(blocker1.
header(
i), chains.param_name(
i));
246 std::stringstream
out;
248 EXPECT_EQ(
"", out.str());
251 ASSERT_EQ(blocker1.
header.size(), chains.num_params());
252 for (
int i = 0;
i < blocker1.
header.size();
i++)
253 EXPECT_EQ(
i, chains.index(blocker1.
header(
i)));
256 std::stringstream
out;
258 EXPECT_EQ(
"", out.str());
262 ASSERT_EQ(1, chains.warmup().size());
263 EXPECT_EQ(0, chains.warmup()(0));
264 EXPECT_EQ(0, chains.warmup(0));
266 chains.set_warmup(10);
267 ASSERT_EQ(1, chains.warmup().size());
268 EXPECT_EQ(10, chains.warmup()(0));
269 EXPECT_EQ(10, chains.warmup(0));
272 chains.set_warmup(100);
273 ASSERT_EQ(1, chains.warmup().size());
274 EXPECT_EQ(100, chains.warmup()(0));
275 EXPECT_EQ(100, chains.warmup(0));
277 chains.add(blocker1);
278 ASSERT_EQ(2, chains.warmup().size());
279 EXPECT_EQ(100, chains.warmup()(0));
280 EXPECT_EQ(100, chains.warmup(0));
281 EXPECT_EQ(0, chains.warmup()(1));
282 EXPECT_EQ(0, chains.warmup(1));
285 chains.set_warmup(1, 20);
286 ASSERT_EQ(2, chains.warmup().size());
287 EXPECT_EQ(100, chains.warmup()(0));
288 EXPECT_EQ(100, chains.warmup(0));
289 EXPECT_EQ(20, chains.warmup()(1));
290 EXPECT_EQ(20, chains.warmup(1));
293 chains.set_warmup(50);
294 ASSERT_EQ(2, chains.warmup().size());
295 EXPECT_EQ(50, chains.warmup()(0));
296 EXPECT_EQ(50, chains.warmup(0));
297 EXPECT_EQ(50, chains.warmup()(1));
298 EXPECT_EQ(50, chains.warmup(1));
301 std::stringstream
out;
304 EXPECT_EQ(
"", out.str());
308 Eigen::VectorXd means1 = blocker1.
samples.colwise().mean();
309 for (
int j = 0;
j < chains.num_params();
j++) {
310 ASSERT_FLOAT_EQ(means1(
j), chains.mean(0,
j))
311 <<
"1: chain, param mean";
312 ASSERT_FLOAT_EQ(means1(
j), chains.mean(
j))
316 chains.add(blocker2);
317 Eigen::VectorXd means2 = blocker2.
samples.colwise().mean();
318 for (
int j = 0;
j < chains.num_params();
j++) {
319 ASSERT_FLOAT_EQ(means2(
j), chains.mean(1,
j))
320 <<
"2: chain, param mean";
321 ASSERT_FLOAT_EQ((means1(
j) + means2(
j)) * 0.5, chains.mean(
j))
326 chains.set_warmup(500);
327 means1 = blocker1.
samples.bottomRows(500).colwise().mean();
328 means2 = blocker2.
samples.bottomRows(500).colwise().mean();
329 for (
int j = 0;
j < chains.num_params();
j++) {
330 ASSERT_FLOAT_EQ(means1(
j), chains.mean(0,
j))
331 <<
"3: chain mean 1 with warmup";
332 ASSERT_FLOAT_EQ(means2(
j), chains.mean(1,
j))
333 <<
"3: chain mean 2 with warmup";
334 ASSERT_FLOAT_EQ((means1(
j) + means2(
j)) * 0.5, chains.mean(
j))
335 <<
"3: param mean with warmup";
338 for (
int j = 0;
j < chains.num_params();
j++) {
340 ASSERT_FLOAT_EQ(chains.mean(0,
j), chains.mean(0,param_name))
341 <<
"4: chain mean 0 called with string name: " << param_name;
342 ASSERT_FLOAT_EQ(chains.mean(1,
j), chains.mean(1,param_name))
343 <<
"4: chain mean 1 called with string name: " << param_name;
344 ASSERT_FLOAT_EQ(chains.mean(
j), chains.mean(param_name))
345 <<
"4: mean called with string name: " << param_name;
349 double sd(Eigen::VectorXd
x) {
350 return sqrt((x.array() - x.mean()).
square().sum() / (x.rows() - 1));
354 std::stringstream
out;
357 EXPECT_EQ(
"", out.str());
362 for (
int j = 0;
j < chains.num_params();
j++) {
363 ASSERT_NEAR(
sd(blocker1.
samples.col(
j)), chains.sd(0,
j), 1
e-8)
364 <<
"1: chain, param sd. index: " <<
j;
365 ASSERT_NEAR(
sd(blocker1.
samples.col(
j)), chains.sd(
j), 1
e-8)
366 <<
"1: param sd. index: " <<
j;
369 chains.add(blocker2);
370 for (
int j = 0;
j < chains.num_params();
j++) {
371 ASSERT_NEAR(
sd(blocker2.
samples.col(
j)), chains.sd(1,
j), 1
e-8)
372 <<
"2: chain, param sd. index: " <<
j;
375 ASSERT_NEAR(
sd(
x), chains.sd(
j), 1
e-8)
376 <<
"2: param sd. index: " <<
j;
379 chains.set_warmup(500);
380 for (
int j = 0;
j < chains.num_params();
j++) {
381 Eigen::VectorXd
x1(500),
x2(500),
x(1000);
382 x1 << blocker1.
samples.col(
j).bottomRows(500);
383 x2 << blocker2.
samples.col(
j).bottomRows(500);
386 ASSERT_NEAR(
sd(x1), chains.sd(0,
j), 1
e-8)
387 <<
"3: chain sd 1 with warmup";
388 ASSERT_NEAR(
sd(x2), chains.sd(1,
j), 1
e-8)
389 <<
"3: chain sd 2 with warmup";
390 ASSERT_NEAR(
sd(x), chains.sd(
j), 1
e-8)
391 <<
"3: param sd with warmup";
394 for (
int j = 0;
j < chains.num_params();
j++) {
396 ASSERT_NEAR(chains.sd(0,
j), chains.sd(0,param_name), 1
e-8)
397 <<
"4: chain sd 0 called with string name: " << param_name;
398 ASSERT_NEAR(chains.sd(1,
j), chains.sd(1,param_name), 1
e-8)
399 <<
"4: chain sd 1 called with string name: " << param_name;
400 ASSERT_NEAR(chains.sd(
j), chains.sd(param_name), 1
e-8)
401 <<
"4: sd called with string name: " << param_name;
408 return (x.array() - x.mean()).
square().sum() / (x.rows() - 1);
412 std::stringstream
out;
415 EXPECT_EQ(
"", out.str());
420 for (
int j = 0;
j < chains.num_params();
j++) {
422 <<
"1: chain, param variance. index: " <<
j;
424 <<
"1: param variance. index: " <<
j;
427 chains.add(blocker2);
428 for (
int j = 0;
j < chains.num_params();
j++) {
430 <<
"2: chain, param variance. index: " <<
j;
434 <<
"2: param variance. index: " <<
j;
437 chains.set_warmup(500);
438 for (
int j = 0;
j < chains.num_params();
j++) {
439 Eigen::VectorXd
x1(500),
x2(500),
x(1000);
440 x1 << blocker1.
samples.col(
j).bottomRows(500);
441 x2 << blocker2.
samples.col(
j).bottomRows(500);
444 ASSERT_NEAR(
variance(x1), chains.variance(0,
j), 1
e-8)
445 <<
"3: chain variance 1 with warmup";
446 ASSERT_NEAR(
variance(x2), chains.variance(1,
j), 1
e-8)
447 <<
"3: chain variance 2 with warmup";
448 ASSERT_NEAR(
variance(x), chains.variance(
j), 1
e-8)
449 <<
"3: param variance with warmup";
452 for (
int j = 0;
j < chains.num_params();
j++) {
454 ASSERT_NEAR(chains.variance(0,
j), chains.variance(0,param_name), 1
e-8)
455 <<
"4: chain variance 0 called with string name: " << param_name;
456 ASSERT_NEAR(chains.variance(1,
j), chains.variance(1,param_name), 1
e-8)
457 <<
"4: chain variance 1 called with string name: " << param_name;
458 ASSERT_NEAR(chains.variance(
j), chains.variance(param_name), 1
e-8)
459 <<
"4: variance called with string name: " << param_name;
465 double x_mean = x.mean();
466 double y_mean = y.mean();
467 return ((x.array() - x_mean) * (y.array() - y_mean)).sum() / (x.rows() - 1);
471 std::stringstream
out;
474 EXPECT_EQ(
"", out.str());
477 chains.add(blocker2);
480 for (
int i = 0;
i < chains.num_params();
i++) {
481 for (
int j =
i;
j < chains.num_params();
j++) {
483 Eigen::VectorXd
x1(1000),
x2(1000),
x(2000);
484 Eigen::VectorXd
y1(1000),
y2(1000),
y(2000);
498 ASSERT_NEAR(cov1, chains.covariance(0,
i,
j), 1
e-8);
499 ASSERT_NEAR(cov2, chains.covariance(1,
i,
j), 1
e-8);
500 ASSERT_NEAR(cov, chains.covariance(
i,
j), 1
e-8);
502 ASSERT_NEAR(cov1, chains.covariance(0,
j,
i), 1
e-8);
503 ASSERT_NEAR(cov2, chains.covariance(1,
j,
i), 1
e-8);
504 ASSERT_NEAR(cov, chains.covariance(
j,
i), 1
e-8);
508 ASSERT_FLOAT_EQ(chains.covariance(0,
i,
j), chains.covariance(0,name1,name2));
509 ASSERT_FLOAT_EQ(chains.covariance(1,
j,
i), chains.covariance(1,name2,name1));
510 ASSERT_FLOAT_EQ(chains.covariance(
i,
j), chains.covariance(name1,name2));
517 std::stringstream
out;
520 EXPECT_EQ(
"", out.str());
523 chains.add(blocker2);
526 for (
int i = 0;
i < chains.num_params();
i++) {
527 for (
int j =
i;
j < chains.num_params();
j++) {
529 Eigen::VectorXd
x1(1000),
x2(1000),
x(2000);
530 Eigen::VectorXd
y1(1000),
y2(1000),
y(2000);
550 corr1 = cov1 /
sd(x1) /
sd(x2);
552 corr2 = cov2 /
sd(y1) /
sd(y2);
554 corr = cov /
sd(x) /
sd(y);
556 ASSERT_NEAR(corr1, chains.correlation(0,
i,
j), 1
e-8)
557 <<
"(" <<
i <<
", " <<
j <<
")";
558 ASSERT_NEAR(corr2, chains.correlation(1,
i,
j), 1
e-8)
559 <<
"(" <<
i <<
", " <<
j <<
")";
560 ASSERT_NEAR(corr, chains.correlation(
i,
j), 1
e-8)
561 <<
"(" <<
i <<
", " <<
j <<
")";
563 ASSERT_NEAR(corr1, chains.correlation(0,
j,
i), 1
e-8)
564 <<
"(" <<
i <<
", " <<
j <<
")";
565 ASSERT_NEAR(corr2, chains.correlation(1,
j,
i), 1
e-8)
566 <<
"(" <<
i <<
", " <<
j <<
")";
567 ASSERT_NEAR(corr, chains.correlation(
j,
i), 1
e-8)
568 <<
"(" <<
i <<
", " <<
j <<
")";
572 ASSERT_FLOAT_EQ(chains.correlation(0,
i,
j), chains.correlation(0,name1,name2));
573 ASSERT_FLOAT_EQ(chains.correlation(1,
j,
i), chains.correlation(1,name2,name1));
574 ASSERT_FLOAT_EQ(chains.correlation(
i,
j), chains.correlation(name1,name2));
581 std::stringstream
out;
583 EXPECT_EQ(
"", out.str());
590 EXPECT_NEAR(0.00241709, chains.quantile(0,index,0.1), 1
e-2);
591 EXPECT_NEAR(0.00348311, chains.quantile(0,index,0.2), 1
e-2);
592 EXPECT_NEAR(0.00477931, chains.quantile(0,index,0.3), 1
e-2);
593 EXPECT_NEAR(0.00607412, chains.quantile(0,index,0.4), 1
e-2);
594 EXPECT_NEAR(0.00770871, chains.quantile(0,index,0.5), 1
e-2);
595 EXPECT_NEAR(0.00999282, chains.quantile(0,index,0.6), 1
e-2);
596 EXPECT_NEAR(0.0131629, chains.quantile(0,index,0.7), 1
e-2);
597 EXPECT_NEAR(0.0185874, chains.quantile(0,index,0.8), 1
e-2);
598 EXPECT_NEAR(0.0263824, chains.quantile(0,index,0.9), 1
e-2);
600 EXPECT_NEAR(0.00241709, chains.quantile(index,0.1), 1
e-2);
601 EXPECT_NEAR(0.00348311, chains.quantile(index,0.2), 1
e-2);
602 EXPECT_NEAR(0.00477931, chains.quantile(index,0.3), 1
e-2);
603 EXPECT_NEAR(0.00607412, chains.quantile(index,0.4), 1
e-2);
604 EXPECT_NEAR(0.00770871, chains.quantile(index,0.5), 1
e-2);
605 EXPECT_NEAR(0.00999282, chains.quantile(index,0.6), 1
e-2);
606 EXPECT_NEAR(0.0131629, chains.quantile(index,0.7), 1
e-2);
607 EXPECT_NEAR(0.0185874, chains.quantile(index,0.8), 1
e-2);
608 EXPECT_NEAR(0.0263824, chains.quantile(index,0.9), 1
e-2);
611 EXPECT_FLOAT_EQ(chains.quantile(0,index,0.1), chains.quantile(0,name,0.1));
612 EXPECT_FLOAT_EQ(chains.quantile(0,index,0.3), chains.quantile(0,name,0.3));
613 EXPECT_FLOAT_EQ(chains.quantile(0,index,0.5), chains.quantile(0,name,0.5));
615 EXPECT_FLOAT_EQ(chains.quantile(index,0.2), chains.quantile(name,0.2));
616 EXPECT_FLOAT_EQ(chains.quantile(index,0.4), chains.quantile(name,0.4));
617 EXPECT_FLOAT_EQ(chains.quantile(index,0.6), chains.quantile(name,0.6));
621 std::stringstream
out;
623 EXPECT_EQ(
"", out.str());
629 Eigen::VectorXd probs(9);
630 probs << 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9;
634 quantiles = chains.quantiles(0,index,probs);
636 ASSERT_EQ(9, quantiles.size());
647 quantiles = chains.quantiles(index,probs);
649 ASSERT_EQ(9, quantiles.size());
662 Eigen::VectorXd quantiles_by_name;
663 quantiles = chains.quantiles(0,index,probs);
664 quantiles_by_name = chains.quantiles(0,name,probs);
666 ASSERT_EQ(quantiles.size(), quantiles_by_name.size());
667 for (
int i = 0;
i < quantiles.size();
i++) {
668 EXPECT_FLOAT_EQ(
quantiles(
i), quantiles_by_name(
i));
671 quantiles = chains.quantiles(index,probs);
672 quantiles_by_name = chains.quantiles(name,probs);
674 ASSERT_EQ(quantiles.size(), quantiles_by_name.size());
675 for (
int i = 0;
i < quantiles.size();
i++) {
676 EXPECT_FLOAT_EQ(
quantiles(
i), quantiles_by_name(
i));
681 std::stringstream
out;
683 EXPECT_EQ(
"", out.str());
691 interval = chains.central_interval(0,index,0.6);
693 EXPECT_NEAR(0.00348311,
interval(0), 1
e-2);
694 EXPECT_NEAR(0.0185874,
interval(1), 1
e-2);
696 interval = chains.central_interval(index,0.6);
698 EXPECT_NEAR(0.00348311,
interval(0), 1
e-2);
699 EXPECT_NEAR(0.0185874,
interval(1), 1
e-2);
702 Eigen::VectorXd interval_by_name;
703 interval = chains.central_interval(0,index,0.6);
704 interval_by_name = chains.central_interval(0,name,0.6);
705 ASSERT_EQ(2, interval_by_name.size());
706 EXPECT_FLOAT_EQ(
interval(0), interval_by_name(0));
707 EXPECT_FLOAT_EQ(
interval(1), interval_by_name(1));
709 interval = chains.central_interval(index,0.6);
710 interval_by_name = chains.central_interval(name,0.6);
711 ASSERT_EQ(2, interval_by_name.size());
712 EXPECT_FLOAT_EQ(
interval(0), interval_by_name(0));
713 EXPECT_FLOAT_EQ(
interval(1), interval_by_name(1));
716 std::stringstream
out;
718 EXPECT_EQ(
"", out.str());
722 EXPECT_NO_THROW(ac = chains.autocorrelation(0,5));
724 EXPECT_NEAR(1, ac[0], 0.01);
725 EXPECT_NEAR(0.529912, ac[1], 0.01);
726 EXPECT_NEAR(0.406604, ac[2], 0.01);
727 EXPECT_NEAR(0.371753, ac[3], 0.01);
728 EXPECT_NEAR(0.310224, ac[4], 0.01);
729 EXPECT_NEAR(0.242701, ac[5], 0.01);
730 EXPECT_NEAR(0.156984, ac[6], 0.01);
731 EXPECT_NEAR(0.112109, ac[7], 0.01);
732 EXPECT_NEAR(0.10186, ac[8], 0.01);
733 EXPECT_NEAR(0.111895, ac[9], 0.01);
734 EXPECT_NEAR(0.117979, ac[10], 0.01);
735 EXPECT_NEAR(0.114381, ac[11], 0.01);
736 EXPECT_NEAR(0.102338, ac[12], 0.01);
737 EXPECT_NEAR(0.108705, ac[13], 0.01);
738 EXPECT_NEAR(0.101822, ac[14], 0.01);
739 EXPECT_NEAR(0.100116, ac[15], 0.01);
740 EXPECT_NEAR(0.110643, ac[16], 0.01);
741 EXPECT_NEAR(0.0732924, ac[17], 0.01);
742 EXPECT_NEAR(0.0500377, ac[18], 0.01);
743 EXPECT_NEAR(0.0221466, ac[19], 0.01);
744 EXPECT_NEAR(0.0548695, ac[20], 0.01);
745 EXPECT_NEAR(0.0778131, ac[21], 0.01);
746 EXPECT_NEAR(0.0618869, ac[22], 0.01);
747 EXPECT_NEAR(0.0734811, ac[23], 0.01);
748 EXPECT_NEAR(0.0719091, ac[24], 0.01);
749 EXPECT_NEAR(0.154124, ac[25], 0.01);
750 EXPECT_NEAR(0.170683, ac[26], 0.01);
751 EXPECT_NEAR(0.0960402, ac[27], 0.01);
752 EXPECT_NEAR(0.140461, ac[28], 0.01);
753 EXPECT_NEAR(0.111866, ac[29], 0.01);
754 EXPECT_NEAR(0.112928, ac[30], 0.01);
757 Eigen::VectorXd ac_by_name;
758 EXPECT_NO_THROW(ac_by_name = chains.autocorrelation(0,name));
759 ASSERT_EQ(ac.size(), ac_by_name.size());
760 for (
int i = 0;
i < ac.size();
i++) {
761 EXPECT_FLOAT_EQ(ac(
i), ac_by_name(
i));
765 std::stringstream
out;
767 EXPECT_EQ(
"", out.str());
771 EXPECT_NO_THROW(ac = chains.autocovariance(0,5));
773 EXPECT_NEAR(0.000150861, ac[0], 0.01);
774 EXPECT_NEAR(7.99431
e-05, ac[1], 0.01);
775 EXPECT_NEAR(6.13408
e-05, ac[2], 0.01);
776 EXPECT_NEAR(5.60831
e-05, ac[3], 0.01);
777 EXPECT_NEAR(4.68008
e-05, ac[4], 0.01);
778 EXPECT_NEAR(3.66142
e-05, ac[5], 0.01);
779 EXPECT_NEAR(2.36828
e-05, ac[6], 0.01);
780 EXPECT_NEAR(1.69129
e-05, ac[7], 0.01);
781 EXPECT_NEAR(1.53667
e-05, ac[8], 0.01);
782 EXPECT_NEAR(1.68806
e-05, ac[9], 0.01);
783 EXPECT_NEAR(1.77985
e-05, ac[10], 0.01);
784 EXPECT_NEAR(1.72556
e-05, ac[11], 0.01);
785 EXPECT_NEAR(1.54389
e-05, ac[12], 0.01);
786 EXPECT_NEAR(1.63994
e-05, ac[13], 0.01);
787 EXPECT_NEAR(1.53609
e-05, ac[14], 0.01);
788 EXPECT_NEAR(1.51037
e-05, ac[15], 0.01);
789 EXPECT_NEAR(1.66917
e-05, ac[16], 0.01);
790 EXPECT_NEAR(1.1057
e-05, ac[17], 0.01);
791 EXPECT_NEAR(7.54875
e-06, ac[18], 0.01);
792 EXPECT_NEAR(3.34107
e-06, ac[19], 0.01);
793 EXPECT_NEAR(8.27767
e-06, ac[20], 0.01);
794 EXPECT_NEAR(1.1739
e-05, ac[21], 0.01);
795 EXPECT_NEAR(9.33633
e-06, ac[22], 0.01);
796 EXPECT_NEAR(1.10854
e-05, ac[23], 0.01);
797 EXPECT_NEAR(1.08483
e-05, ac[24], 0.01);
798 EXPECT_NEAR(2.32514
e-05, ac[25], 0.01);
799 EXPECT_NEAR(2.57494
e-05, ac[26], 0.01);
800 EXPECT_NEAR(1.44887
e-05, ac[27], 0.01);
801 EXPECT_NEAR(2.11901
e-05, ac[28], 0.01);
802 EXPECT_NEAR(1.68763
e-05, ac[29], 0.01);
803 EXPECT_NEAR(1.70365
e-05, ac[30], 0.01);
806 Eigen::VectorXd ac_by_name;
807 EXPECT_NO_THROW(ac_by_name = chains.autocovariance(0,name));
808 ASSERT_EQ(ac.size(), ac_by_name.size());
809 for (
int i = 0;
i < ac.size();
i++) {
810 EXPECT_FLOAT_EQ(ac(
i), ac_by_name(
i));
815 std::stringstream
out;
818 EXPECT_EQ(
"", out.str());
821 chains.add(blocker2);
823 Eigen::VectorXd n_eff(48);
824 n_eff << 466.099,136.953,1170.390,541.256,
825 518.051,589.244,764.813,688.294,
826 323.777,502.892,353.823,588.142,
827 654.336,480.914,176.978,182.649,
828 642.389,470.949,561.947,581.187,
829 446.389,397.641,338.511,678.772,
830 1442.250,837.956,869.865,951.124,
831 619.336,875.805,233.260,786.568,
832 910.144,231.582,907.666,747.347,
833 720.660,195.195,944.547,767.271,
834 723.665,1077.030,470.903,954.924,
835 497.338,583.539,697.204,98.421;
838 ASSERT_NEAR(n_eff(
index - 4), chains.effective_sample_size(
index), 1.0)
839 <<
"n_effective for index: " <<
index <<
", parameter: " 840 << chains.param_name(
index);
845 ASSERT_EQ(chains.effective_sample_size(
index),
846 chains.effective_sample_size(name));
851 std::stringstream
out;
854 EXPECT_EQ(
"", out.str());
857 chains.add(blocker2);
859 Eigen::VectorXd rhat(48);
861 1.00718,1.00473,0.999203,1.00061,1.00378,
862 1.01031,1.00173,1.0045,1.00111,1.00337,
863 1.00546,1.00105,1.00558,1.00463,1.00534,
864 1.01244,1.00174,1.00718,1.00186,1.00554,
865 1.00436,1.00147,1.01017,1.00162,1.00143,
866 1.00058,0.999221,1.00012,1.01028,1.001,
867 1.00305,1.00435,1.00055,1.00246,1.00447,
868 1.0048,1.00209,1.01159,1.00202,1.00077,
869 1.0021,1.00262,1.00308,1.00197,1.00246,
870 1.00085,1.00047,1.00735;
873 ASSERT_NEAR(rhat(
index - 4), chains.split_potential_scale_reduction(
index), 1
e-4)
874 <<
"rhat for index: " <<
index <<
", parameter: " 875 << chains.param_name(
index);
880 ASSERT_EQ(chains.split_potential_scale_reduction(
index),
881 chains.split_potential_scale_reduction(name));
std::vector< double > quantiles(TH1D *h)
const std::string & param_name(int j) const
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
Float_t x1[n_points_granero]
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
double sd(Eigen::VectorXd x)
std::ifstream blocker2_stream
std::ifstream epil2_stream
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
int num_samples(const int chain) const
TEST_F(McmcChains, constructor)
void add(const int chain, const Eigen::MatrixXd &sample)
static stan_csv parse(std::istream &in, std::ostream *out)
double variance(Eigen::VectorXd x)
add("abs", expr_type(int_type()), expr_type(int_type()))
Eigen::Matrix< std::string, Eigen::Dynamic, 1 > header
std::ifstream epil1_stream
const Eigen::VectorXi & warmup() const
std::ifstream blocker1_stream