Dirichlet distribution

dirichlet_distribution.hpp

/*
 * The MIT License (MIT)
 *
 * Copyright (c) 2014 Tsukasa ŌMOTO <henry0312@gmail.com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

/* This file is available under an MIT license. */

#ifndef DIRICHLET_DISTRIBUTION_H
#define DIRICHLET_DISTRIBUTION_H

#include <random>
#include <numeric>
#include <iterator>

/**
 * Dirichlet distribution
 *
 * @see http://en.wikipedia.org/wiki/Dirichlet_distribution
 */
template <class RealType = double>
class dirichlet_distribution
{
public:
    explicit dirichlet_distribution() = default;
    ~dirichlet_distribution() = default;
    template <class Generator, class InputIterator>
        std::vector<RealType> operator()(Generator& gen, const InputIterator first, const InputIterator last);
};

/**
 * Generates samples from Dirichlet distribution parametrized by a vector, alpha (Dir(alpha)).
 *
 * @param first the first iterator of the alpha
 * @param last the last iterator of the alpha
 * @param gen an uniform random number generator object
 * @return a vector (std::vector)
 * @see http://en.wikipedia.org/wiki/Dirichlet_distribution#Gamma_distribution
 */
template <class RealType>
template <class Generator, class InputIterator>
std::vector<RealType> dirichlet_distribution<RealType>::operator()(Generator& gen, const InputIterator first, const InputIterator last)
{
    std::vector<RealType> samples;
    samples.reserve(std::distance(first, last));

    for (InputIterator it = first; it != last; ++it) {
        if (*it > 0) {
            std::gamma_distribution<RealType> dist(*it, 1.0);

            const RealType s = dist(gen);
            if (s < 1.0e-7) {
                // avoid underflow
                samples.push_back(*it); // Mean
            } else {
                samples.push_back(s);
            }
        } else {
            samples.push_back(0.0);
        }
    }

    const RealType sum = std::accumulate(samples.begin(), samples.end(), (RealType)0.0);
    for (auto& s : samples) {
        s /= sum;
    }

    return samples;
}

#endif // DIRICHLET_DISTRIBUTION_H