chains_test.cpp
Go to the documentation of this file.
1 #include <stan/mcmc/chains.hpp>
2 #include <stan/io/stan_csv_reader.hpp>
3 #include <gtest/gtest.h>
4 #include <boost/random/additive_combine.hpp>
5 #include <set>
6 #include <exception>
7 #include <utility>
8 #include <fstream>
9 #include <sstream>
10 #include <vector>
11 #include <string>
12 
13 class McmcChains : public testing::Test {
14 public:
15  void SetUp() {
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");
20  }
21 
22  void TearDown() {
23  blocker1_stream.close();
24  blocker2_stream.close();
25  epil1_stream.close();
26  epil2_stream.close();
27  }
28  std::ifstream blocker1_stream, blocker2_stream, epil1_stream, epil2_stream;
29 };
30 
31 TEST_F(McmcChains, constructor) {
32  std::stringstream out;
33  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
34  EXPECT_EQ("", out.str());
35  // construct with Eigen::Vector
36  stan::mcmc::chains<> chains1(blocker1.header);
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));
41 
42 
43  // construct with stan_csv
44  stan::mcmc::chains<> chains2(blocker1);
45  EXPECT_EQ(1, chains2.num_chains());
46  EXPECT_EQ(blocker1.header.size(), chains2.num_params());
47  for (int i = 0; i < blocker1.header.size(); i++)
48  EXPECT_EQ(blocker1.header(i), chains2.param_name(i));
49  EXPECT_EQ(0, chains2.warmup(0));
50  EXPECT_EQ(1000, chains2.num_samples(0));
51 
52 }
53 
54 TEST_F(McmcChains, add) {
55  std::stringstream out;
56  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
57  EXPECT_EQ("", out.str());
58 
59  // construct with Eigen::Vector
60  stan::mcmc::chains<> chains(blocker1.header);
61  EXPECT_EQ(0, chains.num_chains());
62  EXPECT_EQ(0, chains.num_samples());
63 
64  Eigen::RowVectorXd theta = blocker1.samples.row(0);
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());
71 
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());
79 
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());
88 
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());
97 
98 
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());
108 
109  out.str("");
110  stan::io::stan_csv epil1 = stan::io::stan_csv_reader::parse(epil1_stream, &out);
111  EXPECT_EQ("", out.str());
112  theta.resize(epil1.samples.cols());
113  theta = epil1.samples.row(0);
114  EXPECT_THROW(chains.add(1, theta), std::invalid_argument)
115  << "adding mismatched sample to an existing chain";
116  EXPECT_THROW(chains.add(10, theta), std::invalid_argument)
117  << "adding mismatched sample to a new chain";
118  EXPECT_THROW(chains.add(3, epil1.samples), std::invalid_argument)
119  << "adding mismatched samples to an existing chain";
120  EXPECT_THROW(chains.add(10, epil1.samples), std::invalid_argument)
121  << "adding mismatched samples to a new chain";
122  EXPECT_THROW(chains.add(epil1), std::invalid_argument)
123  << "adding mismatched sample";
124 
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";
139 }
140 
141 TEST_F(McmcChains, add_adapter) {
142  std::stringstream out;
143  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
144  EXPECT_EQ("", out.str());
145 
146  // construct with std::string
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];
150  }
151 
152  stan::mcmc::chains<> chains(param_names);
153  EXPECT_EQ(0, chains.num_chains());
154  EXPECT_EQ(0, chains.num_samples());
155 
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());
159  }
160  for (int i = 0; i < blocker1.samples.rows(); i++) {
161  for (int j = 0; j < blocker1.samples.cols(); j++) {
162  samples[i][j] = blocker1.samples(i, j);
163  }
164  }
165 
166  EXPECT_EQ(blocker1.samples.rows(), static_cast<int>(samples.size()));
167  EXPECT_EQ(blocker1.samples.cols(), static_cast<int>(samples[0].size()));
168 
169  EXPECT_NO_THROW(chains.add(samples))
170  << "adding multiple samples, adds new chain";
171  EXPECT_EQ(1, chains.num_chains());
172  EXPECT_EQ(1000, chains.num_samples(0));
173 
174 }
175 
176 
177 TEST_F(McmcChains, blocker1_num_chains) {
178  std::stringstream out;
179  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
180  EXPECT_EQ("", out.str());
181 
182  stan::mcmc::chains<> chains(blocker1);
183 
184  EXPECT_EQ(1, chains.num_chains());
185 }
186 
187 TEST_F(McmcChains, blocker1_num_samples) {
188  std::stringstream out;
189  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
190  EXPECT_EQ("", out.str());
191 
192  stan::mcmc::chains<> chains(blocker1);
193 
194  EXPECT_EQ(1000, chains.num_samples());
195 }
196 
197 TEST_F(McmcChains, blocker2_num_samples) {
198  std::stringstream out;
199  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
200  EXPECT_EQ("", out.str());
201 
202  stan::mcmc::chains<> chains(blocker2);
203 
204  EXPECT_EQ(1000, chains.num_samples());
205 }
206 
207 TEST_F(McmcChains, blocker_num_samples) {
208  std::stringstream out;
209  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
210  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
211  EXPECT_EQ("", out.str());
212 
213  stan::mcmc::chains<> chains(blocker1);
214  chains.add(blocker2);
215 
216  EXPECT_EQ(2000, chains.num_samples());
217  EXPECT_EQ(1000, chains.num_samples(0));
218  EXPECT_EQ(1000, chains.num_samples(1));
219 }
220 
221 
222 TEST_F(McmcChains, blocker1_param_names) {
223  std::stringstream out;
224  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
225  EXPECT_EQ("", out.str());
226 
227  stan::mcmc::chains<> chains(blocker1);
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));
232  }
233 }
234 TEST_F(McmcChains, blocker1_param_name) {
235  std::stringstream out;
236  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
237  EXPECT_EQ("", out.str());
238 
239  stan::mcmc::chains<> chains(blocker1);
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));
243  }
244 }
245 TEST_F(McmcChains, blocker1_index) {
246  std::stringstream out;
247  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
248  EXPECT_EQ("", out.str());
249 
250  stan::mcmc::chains<> chains(blocker1);
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)));
254 }
255 TEST_F(McmcChains, blocker1_warmup) {
256  std::stringstream out;
257  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
258  EXPECT_EQ("", out.str());
259 
260  stan::mcmc::chains<> chains(blocker1);
261 
262  ASSERT_EQ(1, chains.warmup().size());
263  EXPECT_EQ(0, chains.warmup()(0));
264  EXPECT_EQ(0, chains.warmup(0));
265 
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));
270 
271 
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));
276 
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));
283 
284 
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));
291 
292 
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));
299 }
300 TEST_F(McmcChains, blocker_mean) {
301  std::stringstream out;
302  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
303  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
304  EXPECT_EQ("", out.str());
305 
306  stan::mcmc::chains<> chains(blocker1);
307 
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))
313  << "1: param mean";
314  }
315 
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))
322  << "2: param mean";
323  }
324 
325 
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";
336  }
337 
338  for (int j = 0; j < chains.num_params(); j++) {
339  std::string param_name = chains.param_name(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;
346  }
347 }
348 
349 double sd(Eigen::VectorXd x) {
350  return sqrt((x.array() - x.mean()).square().sum() / (x.rows() - 1));
351 }
352 
353 TEST_F(McmcChains, blocker_sd) {
354  std::stringstream out;
355  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
356  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
357  EXPECT_EQ("", out.str());
358 
359  stan::mcmc::chains<> chains(blocker1);
360 
361  using std::sqrt;
362  for (int j = 0; j < chains.num_params(); j++) {
363  ASSERT_NEAR(sd(blocker1.samples.col(j)), chains.sd(0,j), 1e-8)
364  << "1: chain, param sd. index: " << j;
365  ASSERT_NEAR(sd(blocker1.samples.col(j)), chains.sd(j), 1e-8)
366  << "1: param sd. index: " << j;
367  }
368 
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), 1e-8)
372  << "2: chain, param sd. index: " << j;
373  Eigen::VectorXd x(blocker1.samples.rows() + blocker2.samples.rows());
374  x << blocker1.samples.col(j), blocker2.samples.col(j);
375  ASSERT_NEAR(sd(x), chains.sd(j), 1e-8)
376  << "2: param sd. index: " << j;
377  }
378 
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);
384  x << x1, x2;
385 
386  ASSERT_NEAR(sd(x1), chains.sd(0,j), 1e-8)
387  << "3: chain sd 1 with warmup";
388  ASSERT_NEAR(sd(x2), chains.sd(1,j), 1e-8)
389  << "3: chain sd 2 with warmup";
390  ASSERT_NEAR(sd(x), chains.sd(j), 1e-8)
391  << "3: param sd with warmup";
392  }
393 
394  for (int j = 0; j < chains.num_params(); j++) {
395  std::string param_name = chains.param_name(j);
396  ASSERT_NEAR(chains.sd(0,j), chains.sd(0,param_name), 1e-8)
397  << "4: chain sd 0 called with string name: " << param_name;
398  ASSERT_NEAR(chains.sd(1,j), chains.sd(1,param_name), 1e-8)
399  << "4: chain sd 1 called with string name: " << param_name;
400  ASSERT_NEAR(chains.sd(j), chains.sd(param_name), 1e-8)
401  << "4: sd called with string name: " << param_name;
402  }
403 
404 }
405 
406 
407 double variance(Eigen::VectorXd x) {
408  return (x.array() - x.mean()).square().sum() / (x.rows() - 1);
409 }
410 
411 TEST_F(McmcChains, blocker_variance) {
412  std::stringstream out;
413  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
414  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
415  EXPECT_EQ("", out.str());
416 
417  stan::mcmc::chains<> chains(blocker1);
418 
419  using std::sqrt;
420  for (int j = 0; j < chains.num_params(); j++) {
421  ASSERT_NEAR(variance(blocker1.samples.col(j)), chains.variance(0,j), 1e-8)
422  << "1: chain, param variance. index: " << j;
423  ASSERT_NEAR(variance(blocker1.samples.col(j)), chains.variance(j), 1e-8)
424  << "1: param variance. index: " << j;
425  }
426 
427  chains.add(blocker2);
428  for (int j = 0; j < chains.num_params(); j++) {
429  ASSERT_NEAR(variance(blocker2.samples.col(j)), chains.variance(1,j), 1e-8)
430  << "2: chain, param variance. index: " << j;
431  Eigen::VectorXd x(blocker1.samples.rows() + blocker2.samples.rows());
432  x << blocker1.samples.col(j), blocker2.samples.col(j);
433  ASSERT_NEAR(variance(x), chains.variance(j), 1e-8)
434  << "2: param variance. index: " << j;
435  }
436 
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);
442  x << x1, x2;
443 
444  ASSERT_NEAR(variance(x1), chains.variance(0,j), 1e-8)
445  << "3: chain variance 1 with warmup";
446  ASSERT_NEAR(variance(x2), chains.variance(1,j), 1e-8)
447  << "3: chain variance 2 with warmup";
448  ASSERT_NEAR(variance(x), chains.variance(j), 1e-8)
449  << "3: param variance with warmup";
450  }
451 
452  for (int j = 0; j < chains.num_params(); j++) {
453  std::string param_name = chains.param_name(j);
454  ASSERT_NEAR(chains.variance(0,j), chains.variance(0,param_name), 1e-8)
455  << "4: chain variance 0 called with string name: " << param_name;
456  ASSERT_NEAR(chains.variance(1,j), chains.variance(1,param_name), 1e-8)
457  << "4: chain variance 1 called with string name: " << param_name;
458  ASSERT_NEAR(chains.variance(j), chains.variance(param_name), 1e-8)
459  << "4: variance called with string name: " << param_name;
460  }
461 
462 }
463 
464 double covariance(Eigen::VectorXd x, Eigen::VectorXd y) {
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);
468 }
469 
470 TEST_F(McmcChains, blocker_covariance) {
471  std::stringstream out;
472  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
473  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
474  EXPECT_EQ("", out.str());
475 
476  stan::mcmc::chains<> chains(blocker1);
477  chains.add(blocker2);
478 
479  int n = 0;
480  for (int i = 0; i < chains.num_params(); i++) {
481  for (int j = i; j < chains.num_params(); j++) {
482  if (++n % 13 == 0) { // test every 13th value
483  Eigen::VectorXd x1(1000), x2(1000), x(2000);
484  Eigen::VectorXd y1(1000), y2(1000), y(2000);
485  x1 << blocker1.samples.col(i);
486  x2 << blocker1.samples.col(j);
487 
488  y1 << blocker2.samples.col(i);
489  y2 << blocker2.samples.col(j);
490 
491  x << x1, y1;
492  y << x2, y2;
493 
494  double cov1 = covariance(x1, x2);
495  double cov2 = covariance(y1, y2);
496  double cov = covariance(x, y);
497 
498  ASSERT_NEAR(cov1, chains.covariance(0,i,j), 1e-8);
499  ASSERT_NEAR(cov2, chains.covariance(1,i,j), 1e-8);
500  ASSERT_NEAR(cov, chains.covariance(i,j), 1e-8);
501 
502  ASSERT_NEAR(cov1, chains.covariance(0,j,i), 1e-8);
503  ASSERT_NEAR(cov2, chains.covariance(1,j,i), 1e-8);
504  ASSERT_NEAR(cov, chains.covariance(j,i), 1e-8);
505 
506  std::string name1 = chains.param_name(i);
507  std::string name2 = chains.param_name(j);
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));
511  }
512  }
513  }
514 }
515 
516 TEST_F(McmcChains, blocker_correlation) {
517  std::stringstream out;
518  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
519  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
520  EXPECT_EQ("", out.str());
521 
522  stan::mcmc::chains<> chains(blocker1);
523  chains.add(blocker2);
524 
525  int n = 0;
526  for (int i = 0; i < chains.num_params(); i++) {
527  for (int j = i; j < chains.num_params(); j++) {
528  if (++n % 13 == 0) { // test every 13th value
529  Eigen::VectorXd x1(1000), x2(1000), x(2000);
530  Eigen::VectorXd y1(1000), y2(1000), y(2000);
531  x1 << blocker1.samples.col(i);
532  x2 << blocker1.samples.col(j);
533 
534  y1 << blocker2.samples.col(i);
535  y2 << blocker2.samples.col(j);
536 
537  x << x1, y1;
538  y << x2, y2;
539 
540 
541  double cov1 = covariance(x1, x2);
542  double cov2 = covariance(y1, y2);
543  double cov = covariance(x, y);
544 
545  double corr1 = 0;
546  double corr2 = 0;
547  double corr = 0;
548 
549  if (std::fabs(cov1) > 1e-8)
550  corr1 = cov1 / sd(x1) / sd(x2);
551  if (std::fabs(cov2) > 1e-8)
552  corr2 = cov2 / sd(y1) / sd(y2);
553  if (std::fabs(cov) > 1e-8)
554  corr = cov / sd(x) / sd(y);
555 
556  ASSERT_NEAR(corr1, chains.correlation(0,i,j), 1e-8)
557  << "(" << i << ", " << j << ")";
558  ASSERT_NEAR(corr2, chains.correlation(1,i,j), 1e-8)
559  << "(" << i << ", " << j << ")";
560  ASSERT_NEAR(corr, chains.correlation(i,j), 1e-8)
561  << "(" << i << ", " << j << ")";
562 
563  ASSERT_NEAR(corr1, chains.correlation(0,j,i), 1e-8)
564  << "(" << i << ", " << j << ")";
565  ASSERT_NEAR(corr2, chains.correlation(1,j,i), 1e-8)
566  << "(" << i << ", " << j << ")";
567  ASSERT_NEAR(corr, chains.correlation(j,i), 1e-8)
568  << "(" << i << ", " << j << ")";
569 
570  std::string name1 = chains.param_name(i);
571  std::string name2 = chains.param_name(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));
575  }
576  }
577  }
578 }
579 
580 TEST_F(McmcChains, blocker_quantile) {
581  std::stringstream out;
582  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
583  EXPECT_EQ("", out.str());
584 
585  stan::mcmc::chains<> chains(blocker1);
586 
587  int index = 5;
588 
589  // R's quantile function
590  EXPECT_NEAR(0.00241709, chains.quantile(0,index,0.1), 1e-2);
591  EXPECT_NEAR(0.00348311, chains.quantile(0,index,0.2), 1e-2);
592  EXPECT_NEAR(0.00477931, chains.quantile(0,index,0.3), 1e-2);
593  EXPECT_NEAR(0.00607412, chains.quantile(0,index,0.4), 1e-2);
594  EXPECT_NEAR(0.00770871, chains.quantile(0,index,0.5), 1e-2);
595  EXPECT_NEAR(0.00999282, chains.quantile(0,index,0.6), 1e-2);
596  EXPECT_NEAR(0.0131629, chains.quantile(0,index,0.7), 1e-2);
597  EXPECT_NEAR(0.0185874, chains.quantile(0,index,0.8), 1e-2);
598  EXPECT_NEAR(0.0263824, chains.quantile(0,index,0.9), 1e-2);
599 
600  EXPECT_NEAR(0.00241709, chains.quantile(index,0.1), 1e-2);
601  EXPECT_NEAR(0.00348311, chains.quantile(index,0.2), 1e-2);
602  EXPECT_NEAR(0.00477931, chains.quantile(index,0.3), 1e-2);
603  EXPECT_NEAR(0.00607412, chains.quantile(index,0.4), 1e-2);
604  EXPECT_NEAR(0.00770871, chains.quantile(index,0.5), 1e-2);
605  EXPECT_NEAR(0.00999282, chains.quantile(index,0.6), 1e-2);
606  EXPECT_NEAR(0.0131629, chains.quantile(index,0.7), 1e-2);
607  EXPECT_NEAR(0.0185874, chains.quantile(index,0.8), 1e-2);
608  EXPECT_NEAR(0.0263824, chains.quantile(index,0.9), 1e-2);
609 
610  std::string name = chains.param_name(index);
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));
614 
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));
618 }
619 
620 TEST_F(McmcChains, blocker_quantiles) {
621  std::stringstream out;
622  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
623  EXPECT_EQ("", out.str());
624 
625  stan::mcmc::chains<> chains(blocker1);
626 
627  int index = 5;
628 
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;
631 
632  Eigen::VectorXd quantiles;
633 
634  quantiles = chains.quantiles(0,index,probs);
635  // R's quantile function
636  ASSERT_EQ(9, quantiles.size());
637  EXPECT_NEAR(0.00241709, quantiles(0), 1e-2);
638  EXPECT_NEAR(0.00348311, quantiles(1), 1e-2);
639  EXPECT_NEAR(0.00477931, quantiles(2), 1e-2);
640  EXPECT_NEAR(0.00607412, quantiles(3), 1e-2);
641  EXPECT_NEAR(0.00770871, quantiles(4), 1e-2);
642  EXPECT_NEAR(0.00999282, quantiles(5), 1e-2);
643  EXPECT_NEAR(0.0131629, quantiles(6), 1e-2);
644  EXPECT_NEAR(0.0185874, quantiles(7), 1e-2);
645  EXPECT_NEAR(0.0263824, quantiles(8), 1e-2);
646 
647  quantiles = chains.quantiles(index,probs);
648  // R's quantile function
649  ASSERT_EQ(9, quantiles.size());
650  EXPECT_NEAR(0.00241709, quantiles(0), 1e-2);
651  EXPECT_NEAR(0.00348311, quantiles(1), 1e-2);
652  EXPECT_NEAR(0.00477931, quantiles(2), 1e-2);
653  EXPECT_NEAR(0.00607412, quantiles(3), 1e-2);
654  EXPECT_NEAR(0.00770871, quantiles(4), 1e-2);
655  EXPECT_NEAR(0.00999282, quantiles(5), 1e-2);
656  EXPECT_NEAR(0.0131629, quantiles(6), 1e-2);
657  EXPECT_NEAR(0.0185874, quantiles(7), 1e-2);
658  EXPECT_NEAR(0.0263824, quantiles(8), 1e-2);
659 
660 
661  std::string name = chains.param_name(index);
662  Eigen::VectorXd quantiles_by_name;
663  quantiles = chains.quantiles(0,index,probs);
664  quantiles_by_name = chains.quantiles(0,name,probs);
665 
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));
669  }
670 
671  quantiles = chains.quantiles(index,probs);
672  quantiles_by_name = chains.quantiles(name,probs);
673 
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));
677  }
678 
679 }
680 TEST_F(McmcChains, blocker_central_interval) {
681  std::stringstream out;
682  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
683  EXPECT_EQ("", out.str());
684 
685  stan::mcmc::chains<> chains(blocker1);
686 
687  int index = 5;
688 
689  Eigen::Vector2d interval;
690 
691  interval = chains.central_interval(0,index,0.6);
692  // R's quantile function
693  EXPECT_NEAR(0.00348311, interval(0), 1e-2); // 0.2
694  EXPECT_NEAR(0.0185874, interval(1), 1e-2); // 0.8
695 
696  interval = chains.central_interval(index,0.6);
697  // R's quantile function
698  EXPECT_NEAR(0.00348311, interval(0), 1e-2); // 0.2
699  EXPECT_NEAR(0.0185874, interval(1), 1e-2); // 0.8
700 
701  std::string name = chains.param_name(index);
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));
708 
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));
714 }
715 TEST_F(McmcChains, blocker_autocorrelation) {
716  std::stringstream out;
717  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
718  EXPECT_EQ("", out.str());
719 
720  stan::mcmc::chains<> chains(blocker1);
721  Eigen::VectorXd ac;
722  EXPECT_NO_THROW(ac = chains.autocorrelation(0,5));
723 
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);
755 
756  std::string name = chains.param_name(5);
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));
762  }
763 }
764 TEST_F(McmcChains, blocker_autocovariance) {
765  std::stringstream out;
766  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
767  EXPECT_EQ("", out.str());
768 
769  stan::mcmc::chains<> chains(blocker1);
770  Eigen::VectorXd ac;
771  EXPECT_NO_THROW(ac = chains.autocovariance(0,5));
772 
773  EXPECT_NEAR(0.000150861, ac[0], 0.01);
774  EXPECT_NEAR(7.99431e-05, ac[1], 0.01);
775  EXPECT_NEAR(6.13408e-05, ac[2], 0.01);
776  EXPECT_NEAR(5.60831e-05, ac[3], 0.01);
777  EXPECT_NEAR(4.68008e-05, ac[4], 0.01);
778  EXPECT_NEAR(3.66142e-05, ac[5], 0.01);
779  EXPECT_NEAR(2.36828e-05, ac[6], 0.01);
780  EXPECT_NEAR(1.69129e-05, ac[7], 0.01);
781  EXPECT_NEAR(1.53667e-05, ac[8], 0.01);
782  EXPECT_NEAR(1.68806e-05, ac[9], 0.01);
783  EXPECT_NEAR(1.77985e-05, ac[10], 0.01);
784  EXPECT_NEAR(1.72556e-05, ac[11], 0.01);
785  EXPECT_NEAR(1.54389e-05, ac[12], 0.01);
786  EXPECT_NEAR(1.63994e-05, ac[13], 0.01);
787  EXPECT_NEAR(1.53609e-05, ac[14], 0.01);
788  EXPECT_NEAR(1.51037e-05, ac[15], 0.01);
789  EXPECT_NEAR(1.66917e-05, ac[16], 0.01);
790  EXPECT_NEAR(1.1057e-05, ac[17], 0.01);
791  EXPECT_NEAR(7.54875e-06, ac[18], 0.01);
792  EXPECT_NEAR(3.34107e-06, ac[19], 0.01);
793  EXPECT_NEAR(8.27767e-06, ac[20], 0.01);
794  EXPECT_NEAR(1.1739e-05, ac[21], 0.01);
795  EXPECT_NEAR(9.33633e-06, ac[22], 0.01);
796  EXPECT_NEAR(1.10854e-05, ac[23], 0.01);
797  EXPECT_NEAR(1.08483e-05, ac[24], 0.01);
798  EXPECT_NEAR(2.32514e-05, ac[25], 0.01);
799  EXPECT_NEAR(2.57494e-05, ac[26], 0.01);
800  EXPECT_NEAR(1.44887e-05, ac[27], 0.01);
801  EXPECT_NEAR(2.11901e-05, ac[28], 0.01);
802  EXPECT_NEAR(1.68763e-05, ac[29], 0.01);
803  EXPECT_NEAR(1.70365e-05, ac[30], 0.01);
804 
805  std::string name = chains.param_name(5);
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));
811  }
812 }
813 
814 TEST_F(McmcChains,blocker_effective_sample_size) {
815  std::stringstream out;
816  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
817  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
818  EXPECT_EQ("", out.str());
819 
820  stan::mcmc::chains<> chains(blocker1);
821  chains.add(blocker2);
822 
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;
836 
837  for (int index = 4; index < chains.num_params(); index++) {
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);
841  }
842 
843  for (int index = 0; index < chains.num_params(); index++) {
844  std::string name = chains.param_name(index);
845  ASSERT_EQ(chains.effective_sample_size(index),
846  chains.effective_sample_size(name));
847  }
848 }
849 
850 TEST_F(McmcChains,blocker_split_potential_scale_reduction) {
851  std::stringstream out;
852  stan::io::stan_csv blocker1 = stan::io::stan_csv_reader::parse(blocker1_stream, &out);
853  stan::io::stan_csv blocker2 = stan::io::stan_csv_reader::parse(blocker2_stream, &out);
854  EXPECT_EQ("", out.str());
855 
856  stan::mcmc::chains<> chains(blocker1);
857  chains.add(blocker2);
858 
859  Eigen::VectorXd rhat(48);
860  rhat <<
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;
871 
872  for (int index = 4; index < chains.num_params(); index++) {
873  ASSERT_NEAR(rhat(index - 4), chains.split_potential_scale_reduction(index), 1e-4)
874  << "rhat for index: " << index << ", parameter: "
875  << chains.param_name(index);
876  }
877 
878  for (int index = 0; index < chains.num_params(); index++) {
879  std::string name = chains.param_name(index);
880  ASSERT_EQ(chains.split_potential_scale_reduction(index),
881  chains.split_potential_scale_reduction(name));
882  }
883 
884 }
name
const XML_Char * name
Definition: expat.h:151
quantiles
std::vector< double > quantiles(TH1D *h)
Definition: absCal.cxx:528
stan::mcmc::chains::param_name
const std::string & param_name(int j) const
Definition: chains.hpp:368
stan::io::stan_csv::samples
Eigen::MatrixXd samples
Definition: stan_csv_reader.hpp:65
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
McmcChains::TearDown
void TearDown()
Definition: chains_test.cpp:22
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
submit_syst.x2
x2
Definition: submit_syst.py:250
covariance
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
Definition: chains_test.cpp:464
submit_syst.y2
y2
Definition: submit_syst.py:250
calib::corr
double corr(TGraph *g, double thresh)
Definition: TestParticleCorrections_module.cc:110
stan_csv_reader.hpp
cet::square
constexpr T square(T x)
Definition: pow.h:21
submit_syst.y
y
Definition: submit_syst.py:249
McmcChains::SetUp
void SetUp()
Definition: chains_test.cpp:15
sd
double sd(Eigen::VectorXd x)
Definition: chains_test.cpp:349
submit_syst.x
x
Definition: submit_syst.py:249
McmcChains::blocker2_stream
std::ifstream blocker2_stream
Definition: chains_test.cpp:28
McmcChains::epil2_stream
std::ifstream epil2_stream
Definition: chains_test.cpp:28
stan::io::stan_csv
Definition: stan_csv_reader.hpp:61
calib::j
const double j
Definition: BetheBloch.cxx:29
stan::math::invalid_argument
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Definition: invalid_argument.hpp:30
stan::mcmc::chains::num_samples
int num_samples(const int chain) const
Definition: chains.hpp:396
TEST_F
TEST_F(McmcChains, constructor)
Definition: chains_test.cpp:31
xml_schema::string
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
stan::mcmc::chains::num_params
int num_params() const
Definition: chains.hpp:360
stan::mcmc::chains::add
void add(const int chain, const Eigen::MatrixXd &sample)
Definition: chains.hpp:418
stan::mcmc::chains::num_chains
int num_chains() const
Definition: chains.hpp:356
stan::io::stan_csv_reader::parse
static stan_csv parse(std::istream &in, std::ostream *out)
Definition: stan_csv_reader.hpp:329
e
Float_t e
Definition: plot.C:35
variance
double variance(Eigen::VectorXd x)
Definition: chains_test.cpp:407
add
add("abs", expr_type(int_type()), expr_type(int_type()))
stan::io::stan_csv::header
Eigen::Matrix< std::string, Eigen::Dynamic, 1 > header
Definition: stan_csv_reader.hpp:63
McmcChains::epil1_stream
std::ifstream epil1_stream
Definition: chains_test.cpp:28
stan::mcmc::chains::warmup
const Eigen::VectorXi & warmup() const
Definition: chains.hpp:388
McmcChains::blocker1_stream
std::ifstream blocker1_stream
Definition: chains_test.cpp:28
Generated on Sat Sep 26 2020 05:14:55 for NOvA by   doxygen 1.8.11