Simulation of simple discrete distributions

distributions.cpp

#include <chrono>
#include <iostream>
#include <random>
#include <vector>

using std::cout;
using std::endl;
using std::vector;

class TestRandom {
public:
        TestRandom() {
        seed = std::chrono::system_clock::now().time_since_epoch().count();
        generator.seed(seed);
    }

    bool bernoulliDistribution(double p);
    void testBernoulliDistribution();

    int binomialDistribution(int n, double p);
    void testBinomialDistribution();

    int geometricDistribution(double p);
    void testGeometricDistribution();

private:
    int seed;
    std::minstd_rand generator;
    const int sampleCount = 10000;

    double meanValue(const vector<int>& v) {
        auto meanValue = 0.0;
        for(auto e : v)
            meanValue += e;
        meanValue /= v.size();
        return meanValue;
    }

    double variance(const vector<int>& v) {
        auto variance = 0.0;
        auto m = meanValue(v);
        for (auto e : v)
            variance += pow(e - m, 2);
        variance /= v.size();
        return variance;
    }

};

bool TestRandom::bernoulliDistribution(double p) {
    return generator() < generator.min() + p*(generator.max() - generator.min());
}

void TestRandom::testBernoulliDistribution() {
    auto p = 0.3;
    cout << "----------------------------------------" << endl;
    cout << "Test for Bernoulli distribution with p=" << p << endl;
    auto failCount = 0;
    auto successCount = 0;
    vector<int> sample(sampleCount, 0);
    for (auto& e : sample) {
        e = bernoulliDistribution(p);
        e ? successCount++ : failCount++;
    }
    cout << "Success count= " << successCount << endl;
    cout << "Fail count= " << failCount << endl;
    cout << "Success rate= " << 
        static_cast<double>(successCount) / sampleCount << endl;
    cout << "Fail rate= " << 
        static_cast<double>(failCount) / sampleCount << endl;
    cout << "Sample mean value= " << meanValue(sample) << endl;
    cout << "Sample variance= " << variance(sample) << endl;
    cout << "Theoretical mean value: p=" << p << endl;
    cout << "Theoretical variance: p*(1-p)=" << p*(1-p) << endl;
}

int TestRandom::binomialDistribution(int n, double p) {
    int successRate = 0;
    int failRate = 0;
    for (int i = 0; i < n; i++) {
        bernoulliDistribution(p) ? successRate++ : failRate++;
    }
    return successRate;
}

void TestRandom::testBinomialDistribution() {
    auto n = 10;
    auto p = 0.7;
    cout << "----------------------------------------" << endl;
    cout << "Testing of Binomial distribution(" << n <<
        "," << p << ")" << endl;

    vector<int> sample(sampleCount, 0);
    for(auto& e : sample)
        e = binomialDistribution(n, p);

    cout << "Resulting hystogram" << endl;
    vector<double> hyst(n + 1, 0);
    for (auto e : sample) {
        hyst[e]++;
    }

    for (auto& e : hyst) {
        cout << (static_cast<double>(e)*100 / sampleCount) << " ";
    }
    cout << endl;

    auto sampleMeanValue = meanValue(sample);
    auto sampleVariance = variance(sample);

    cout << "Sample mean value: " << sampleMeanValue << endl;
    cout << "Sample variance: " << sampleVariance << endl;

    cout << "Theoretical mean value: n*p=" << n*p << endl;
    cout << "Theoretical variance:  n*p*(1-p)=" << n*p*(1 - p) << endl;
}

int TestRandom::geometricDistribution(double p) {
    auto r = 0;
    while (bernoulliDistribution(p) == 0)
        r++;
    return r + 1;
}

void TestRandom::testGeometricDistribution() {
    auto p = 0.3;
    cout << "----------------------------------------" << endl;
    cout << "Testing of Geometric distribution(" << p << ")" << endl;
    vector<int> sample(sampleCount, 0);
    for (auto& e : sample) {
        e = geometricDistribution(p);
    }
    auto sampleMeanValue = meanValue(sample);
    auto sampleVariance = variance(sample);

    cout << "Sample mean value: " << sampleMeanValue << endl;
    cout << "Sample variance: " << sampleVariance << endl;

    cout << "Theoretical mean value: 1/p=" << 1.0/p << endl;
    cout << "Theoretical variance:  (1-p)/p^2=" << (1-p)/(p*p) << endl;
}

int main()
{
    TestRandom tr;
    tr.testBernoulliDistribution();
    tr.testBinomialDistribution();
    tr.testGeometricDistribution();
    return 0;
}