gistfile1.cpp

#include <algorithm>
#include <iostream>
#include <iterator>
#include <random>
#include <vector>
using namespace std;

class multinominal_distribution
{
    std::vector<double> m_prob;
    std::mt19937_64     m_eng;
    std::uniform_real_distribution<double>
                        m_dist;

public:
    template <typename T>
    multinominal_distribution(std::initializer_list<T> prob)
        : m_prob(prob.begin(), prob.end())
        , m_dist(0, std::accumulate(prob.begin(), prob.end(), 0.0))
    {
    }

    size_t operator()()
    {
        double rand = m_dist(m_eng);
        double sum = 0;
        for (size_t k = 0; k < m_prob.size(); ++k) {
            sum += m_prob[k];
            if (rand < sum) return k;
        }
        return m_prob.size() - 1;
    }
};

int main()
{
    multinominal_distribution multi({ 1, 3, 6 }); 

    vector<int> acc(3, 0);
    for (int i = 0; i < 100000; ++i) {
        acc[multi()] += 1;
    }

    copy(acc.begin(), acc.end(), ostream_iterator<int>(cout,","));
}

gistfile2.cpp

#include <iostream>
#include <iterator>
#include <random>
using namespace std;

template <typename Rng, typename T> double accumulate(Rng const& rng, T init)
{
    return std::accumulate(std::begin(rng), std::end(rng), init);
}

class dirichret_distribution {
    std::vector<std::mt19937_64> m_eng;
    std::vector<std::gamma_distribution<double> > m_dist;

  public:
    template <typename T> dirichret_distribution(std::initializer_list<T> alpha)
    {
        for (auto& a : alpha) {
            m_eng.emplace_back();
            m_dist.emplace_back(a, 1.0);
        }
    }

    std::vector<double> operator()()
    {
        std::vector<double> rand;
        for (size_t i = 0; i < m_dist.size(); ++i) {
            rand.push_back(m_dist[i](m_eng[i]));
        }
        double sum = accumulate(rand, 0.0);
        if (sum > 0)
            for (auto& r : rand)
                r /= sum;
        return rand;
    }
};

int main(int argc, char const* argv[])
{
    dirichret_distribution dir({ 1.0, 0.1, 0.1 });

    vector<vector<double> > acc;
    for (size_t i = 0; i < 1000; ++i) {
        acc.push_back(dir());
    }

    vector<double> mean;
    for (int i = 0; i < 3; ++i) {
        double sum = 0;
        for (int k = 0; k < acc.size(); ++k)
            sum += acc[k][i];
        mean.push_back(sum / acc.size());
    }

    copy(mean.begin(), mean.end(), ostream_iterator<double>(cout, ","));
}