1 /* boost histogram.cpp graphical verification of distribution functions
2 *
3 * Copyright Jens Maurer 2000
4 * Distributed under the Boost Software License, Version 1.0. (See
5 * accompanying file LICENSE_1_0.txt or copy at
6 * http://www.boost.org/LICENSE_1_0.txt)
7 *
8 * $Id$
9 *
10 * This test program allows to visibly examine the results of the
11 * distribution functions.
12 */
13
14 #include <iostream>
15 #include <iomanip>
16 #include <vector>
17 #include <algorithm>
18 #include <cmath>
19 #include <string>
20 #include <boost/random.hpp>
21
22
plot_histogram(const std::vector<int> & slots,int samples,double from,double to)23 void plot_histogram(const std::vector<int>& slots, int samples,
24 double from, double to)
25 {
26 int m = *std::max_element(slots.begin(), slots.end());
27 const int nRows = 20;
28 std::cout.setf(std::ios::fixed|std::ios::left);
29 std::cout.precision(5);
30 for(int r = 0; r < nRows; r++) {
31 double y = ((nRows - r) * double(m))/(nRows * samples);
32 std::cout << std::setw(10) << y << " ";
33 for(unsigned int col = 0; col < slots.size(); col++) {
34 char out = ' ';
35 if(slots[col]/double(samples) >= y)
36 out = 'x';
37 std::cout << out;
38 }
39 std::cout << std::endl;
40 }
41 std::cout << std::setw(12) << " "
42 << std::setw(10) << from;
43 std::cout.setf(std::ios::right, std::ios::adjustfield);
44 std::cout << std::setw(slots.size()-10) << to << std::endl;
45 }
46
47 // I am not sure whether these two should be in the library as well
48
49 // maintain sum of NumberGenerator results
50 template<class NumberGenerator,
51 class Sum = typename NumberGenerator::result_type>
52 class sum_result
53 {
54 public:
55 typedef NumberGenerator base_type;
56 typedef typename base_type::result_type result_type;
sum_result(const base_type & g)57 explicit sum_result(const base_type & g) : gen(g), _sum(0) { }
operator ()()58 result_type operator()() { result_type r = gen(); _sum += r; return r; }
base()59 base_type & base() { return gen; }
sum() const60 Sum sum() const { return _sum; }
reset()61 void reset() { _sum = 0; }
62 private:
63 base_type gen;
64 Sum _sum;
65 };
66
67
68 // maintain square sum of NumberGenerator results
69 template<class NumberGenerator,
70 class Sum = typename NumberGenerator::result_type>
71 class squaresum_result
72 {
73 public:
74 typedef NumberGenerator base_type;
75 typedef typename base_type::result_type result_type;
squaresum_result(const base_type & g)76 explicit squaresum_result(const base_type & g) : gen(g), _sum(0) { }
operator ()()77 result_type operator()() { result_type r = gen(); _sum += r*r; return r; }
base()78 base_type & base() { return gen; }
squaresum() const79 Sum squaresum() const { return _sum; }
reset()80 void reset() { _sum = 0; }
81 private:
82 base_type gen;
83 Sum _sum;
84 };
85
86
87 template<class RNG>
histogram(RNG base,int samples,double from,double to,const std::string & name)88 void histogram(RNG base, int samples, double from, double to,
89 const std::string & name)
90 {
91 typedef squaresum_result<sum_result<RNG, double>, double > SRNG;
92 SRNG gen((sum_result<RNG, double>(base)));
93 const int nSlots = 60;
94 std::vector<int> slots(nSlots,0);
95 for(int i = 0; i < samples; i++) {
96 double val = gen();
97 if(val < from || val >= to) // early check avoids overflow
98 continue;
99 int slot = int((val-from)/(to-from) * nSlots);
100 if(slot < 0 || slot > (int)slots.size())
101 continue;
102 slots[slot]++;
103 }
104 std::cout << name << std::endl;
105 plot_histogram(slots, samples, from, to);
106 double mean = gen.base().sum() / samples;
107 std::cout << "mean: " << mean
108 << " sigma: " << std::sqrt(gen.squaresum()/samples-mean*mean)
109 << "\n" << std::endl;
110 }
111
112 template<class PRNG, class Dist>
make_gen(PRNG & rng,Dist d)113 inline boost::variate_generator<PRNG&, Dist> make_gen(PRNG & rng, Dist d)
114 {
115 return boost::variate_generator<PRNG&, Dist>(rng, d);
116 }
117
118 template<class PRNG>
histograms()119 void histograms()
120 {
121 PRNG rng;
122 using namespace boost;
123 histogram(make_gen(rng, uniform_smallint<>(0, 5)), 100000, -1, 6,
124 "uniform_smallint(0,5)");
125 histogram(make_gen(rng, uniform_int<>(0, 5)), 100000, -1, 6,
126 "uniform_int(0,5)");
127 histogram(make_gen(rng, uniform_real<>(0,1)), 100000, -0.5, 1.5,
128 "uniform_real(0,1)");
129 histogram(make_gen(rng, bernoulli_distribution<>(0.2)), 100000, -0.5, 1.5,
130 "bernoulli(0.2)");
131 histogram(make_gen(rng, binomial_distribution<>(4, 0.2)), 100000, -1, 5,
132 "binomial(4, 0.2)");
133 histogram(make_gen(rng, triangle_distribution<>(1, 2, 8)), 100000, 0, 10,
134 "triangle(1,2,8)");
135 histogram(make_gen(rng, geometric_distribution<>(5.0/6.0)), 100000, 0, 10,
136 "geometric(5/6)");
137 histogram(make_gen(rng, exponential_distribution<>(0.3)), 100000, 0, 10,
138 "exponential(0.3)");
139 histogram(make_gen(rng, cauchy_distribution<>()), 100000, -5, 5,
140 "cauchy");
141 histogram(make_gen(rng, lognormal_distribution<>(3, 2)), 100000, 0, 10,
142 "lognormal");
143 histogram(make_gen(rng, normal_distribution<>()), 100000, -3, 3,
144 "normal");
145 histogram(make_gen(rng, normal_distribution<>(0.5, 0.5)), 100000, -3, 3,
146 "normal(0.5, 0.5)");
147 histogram(make_gen(rng, poisson_distribution<>(1.5)), 100000, 0, 5,
148 "poisson(1.5)");
149 histogram(make_gen(rng, poisson_distribution<>(10)), 100000, 0, 20,
150 "poisson(10)");
151 histogram(make_gen(rng, gamma_distribution<>(0.5)), 100000, 0, 0.5,
152 "gamma(0.5)");
153 histogram(make_gen(rng, gamma_distribution<>(1)), 100000, 0, 3,
154 "gamma(1)");
155 histogram(make_gen(rng, gamma_distribution<>(2)), 100000, 0, 6,
156 "gamma(2)");
157 }
158
159
main()160 int main()
161 {
162 histograms<boost::mt19937>();
163 // histograms<boost::lagged_fibonacci607>();
164 }
165
166