1 | #include <boost/random.hpp>
|
---|
2 | #include <boost/math/distributions.hpp>
|
---|
3 | #include <vector>
|
---|
4 | #include <utility>
|
---|
5 | #include <math.h>
|
---|
6 | #include <iostream>
|
---|
7 |
|
---|
8 | std::vector<int> create_samples(const boost::random::negative_binomial_distribution<> &nb, const int size)
|
---|
9 | {
|
---|
10 | boost::random::mt19937 rng;
|
---|
11 |
|
---|
12 | std::vector<int> v(size);
|
---|
13 | for(std::vector<int>::iterator it = v.begin(); it != v.end(); ++it)
|
---|
14 | *it = nb(rng);
|
---|
15 |
|
---|
16 | return v;
|
---|
17 | }
|
---|
18 |
|
---|
19 | std::pair<double,double> mean_stddev(const std::vector<int> &v)
|
---|
20 | {
|
---|
21 | double n_elements = static_cast<double>(v.size());
|
---|
22 | double sum = std::accumulate(std::begin(v), std::end(v), 0.0);
|
---|
23 | double m = sum / n_elements;
|
---|
24 |
|
---|
25 | double accum = 0.0;
|
---|
26 | std::for_each (std::begin(v), std::end(v), [&](const double d) {
|
---|
27 | accum += (d - m) * (d - m);
|
---|
28 | });
|
---|
29 |
|
---|
30 | double stdev = sqrt(accum / (n_elements-1));
|
---|
31 |
|
---|
32 | return std::make_pair(m,stdev);
|
---|
33 | }
|
---|
34 |
|
---|
35 | int main()
|
---|
36 | {
|
---|
37 | const int size = 1000000;
|
---|
38 | double _k = 2.99;
|
---|
39 | double _p = 0.01;
|
---|
40 |
|
---|
41 | boost::random::negative_binomial_distribution<> nb_sample(_k,_p);
|
---|
42 |
|
---|
43 | std::vector<int> v = create_samples(nb_sample, size);
|
---|
44 | std::vector<int>::const_iterator it = v.begin();
|
---|
45 |
|
---|
46 |
|
---|
47 | std::pair<double,double> ms = mean_stddev(v);
|
---|
48 | std::cout << std::endl << "after " << size << " simulations of negative binomial distribution (k=" << _k << ", p=" << _p << "): " << std::endl;
|
---|
49 | std::cout << " mean = " << ms.first << std::endl;
|
---|
50 | std::cout << "stdev = " << ms.second << std::endl;
|
---|
51 |
|
---|
52 | boost::math::negative_binomial_distribution<> nb(_k,_p);
|
---|
53 | std::cout << std::endl << "boost::math::negative_binomial nb(" << nb.successes() << "," << nb.success_fraction() << ") gives:" << std::endl;
|
---|
54 | std::cout << " mean = " << boost::math::mean(nb) << std::endl;
|
---|
55 | std::cout << "stdev = " << std::sqrt(boost::math::variance(nb)) << std::endl;
|
---|
56 |
|
---|
57 | boost::math::negative_binomial_distribution<> nb_((int)_k,_p);
|
---|
58 | std::cout << std::endl << "boost::math::negative_binomial nb(" << nb_.successes() << "," << nb_.success_fraction() << ") gives:" << std::endl;
|
---|
59 | std::cout << " mean = " << boost::math::mean(nb_) << std::endl;
|
---|
60 | std::cout << "stdev = " << std::sqrt(boost::math::variance(nb_)) << std::endl;
|
---|
61 |
|
---|
62 | return 0;
|
---|
63 | }
|
---|