[29] | 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: histogram.cpp,v 1.8 2004/07/27 03:43:34 dgregor Exp $ |
---|
| 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 | |
---|
| 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; |
---|
| 57 | explicit sum_result(const base_type & g) : gen(g), _sum(0) { } |
---|
| 58 | result_type operator()() { result_type r = gen(); _sum += r; return r; } |
---|
| 59 | base_type & base() { return gen; } |
---|
| 60 | Sum sum() const { return _sum; } |
---|
| 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; |
---|
| 76 | explicit squaresum_result(const base_type & g) : gen(g), _sum(0) { } |
---|
| 77 | result_type operator()() { result_type r = gen(); _sum += r*r; return r; } |
---|
| 78 | base_type & base() { return gen; } |
---|
| 79 | Sum squaresum() const { return _sum; } |
---|
| 80 | void reset() { _sum = 0; } |
---|
| 81 | private: |
---|
| 82 | base_type gen; |
---|
| 83 | Sum _sum; |
---|
| 84 | }; |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | template<class RNG> |
---|
| 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> |
---|
| 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> |
---|
| 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 | |
---|
| 160 | int main() |
---|
| 161 | { |
---|
| 162 | histograms<boost::mt19937>(); |
---|
| 163 | // histograms<boost::lagged_fibonacci607>(); |
---|
| 164 | } |
---|
| 165 | |
---|