#include <boost/random.hpp>
#include <boost/math/distributions.hpp>
#include <vector>
#include <utility>
#include <math.h> 
#include <iostream>

std::vector<int> create_samples(const boost::random::negative_binomial_distribution<> &nb, const int size)
{
    boost::random::mt19937 rng;

    std::vector<int> v(size);
    for(std::vector<int>::iterator it = v.begin(); it != v.end(); ++it)
        *it = nb(rng);

    return v;
}

std::pair<double,double> mean_stddev(const std::vector<int> &v)
{
    double n_elements = static_cast<double>(v.size());
    double sum = std::accumulate(std::begin(v), std::end(v), 0.0);
    double m =  sum / n_elements;

    double accum = 0.0;
    std::for_each (std::begin(v), std::end(v), [&](const double d) {
        accum += (d - m) * (d - m);
    });

    double stdev = sqrt(accum / (n_elements-1));

    return std::make_pair(m,stdev);
}

int main()
{
    const int size = 1000000;
    double _k = 2.99;
    double _p = 0.01;

    boost::random::negative_binomial_distribution<> nb_sample(_k,_p);

    std::vector<int> v = create_samples(nb_sample, size);
    std::vector<int>::const_iterator it = v.begin();


    std::pair<double,double> ms = mean_stddev(v);
    std::cout << std::endl << "after " << size << " simulations of negative binomial distribution (k=" << _k << ", p=" << _p << "): " << std::endl;
    std::cout << " mean = " << ms.first << std::endl;
    std::cout << "stdev = " << ms.second << std::endl;

    boost::math::negative_binomial_distribution<> nb(_k,_p);
    std::cout << std::endl << "boost::math::negative_binomial nb(" << nb.successes() << "," << nb.success_fraction() << ") gives:" << std::endl;
    std::cout << " mean = " << boost::math::mean(nb) << std::endl;
    std::cout << "stdev = " << std::sqrt(boost::math::variance(nb)) << std::endl;

    boost::math::negative_binomial_distribution<> nb_((int)_k,_p);
    std::cout << std::endl << "boost::math::negative_binomial nb(" << nb_.successes() << "," << nb_.success_fraction() << ") gives:" << std::endl;
    std::cout << " mean = " << boost::math::mean(nb_) << std::endl;
    std::cout << "stdev = " << std::sqrt(boost::math::variance(nb_)) << std::endl;

    return 0;
}