Probability.hpp

Go to the documentation of this file.

#ifndef AI_TOOLBOX_UTILS_PROBABILITY_HEADER_FILE
#define AI_TOOLBOX_UTILS_PROBABILITY_HEADER_FILE
 
#include <random>
#include <algorithm>
 
#include <AIToolbox/Types.hpp>
#include <AIToolbox/Utils/Core.hpp>
 
namespace AIToolbox {
    static std::uniform_real_distribution<double> probabilityDistribution(0.0, 1.0);
 
    template <typename T>
    bool isProbability(const size_t size, const T & in) {
        double p = 0.0;
        for (size_t i = 0; i < size; ++i) {
            const double value = static_cast<double>(in[i]);
            if (value < 0.0) return false;
            p += value;
        }
        if (checkDifferentSmall(p, 1.0))
            return false;
 
        return true;
    }
 
    template <typename T>
    bool isProbability(const size_t rows, const size_t cols, const T & in) {
        for (size_t row = 0; row < rows; ++row)
            if (!isProbability(cols, in[row]))
                return false;
        return true;
    }
 
    template <typename T>
    bool isProbability(const size_t depth, const size_t rows, const size_t cols, const T & in) {
        for (size_t d = 0; d < depth; ++d)
            if (!isProbability(rows, cols, in[d]))
                return false;
        return true;
    }
 
    bool isProbability(const Matrix2D & in);
 
    bool isProbability(const Matrix3D & in);
 
    bool isProbability(const SparseMatrix2D & in);
 
    bool isProbability(const SparseMatrix3D & in);
 
    template <typename T, typename G>
    size_t sampleProbability(const size_t d, const T& in, G& generator) {
        double p = probabilityDistribution(generator);
 
        for ( size_t i = 0; i < d; ++i ) {
            if ( in[i] > p ) return i;
            p -= in[i];
        }
        return d-1;
    }
 
    template <typename G>
    size_t sampleProbability(const size_t d, const SparseMatrix2D::ConstRowXpr& in, G& generator) {
        double p = probabilityDistribution(generator);
 
        for ( SparseMatrix2D::ConstRowXpr::InnerIterator i(in, 0); ; ++i ) {
            if ( i.value() > p ) return i.col();
            p -= i.value();
        }
        return d-1;
    }
 
    template <typename G>
    double sampleBetaDistribution(double a, double b, G & generator) {
        std::gamma_distribution<double> dista(a, 1.0);
        std::gamma_distribution<double> distb(b, 1.0);
        const auto X = dista(generator);
        const auto Y = distb(generator);
        return X / (X + Y);
    }
 
    template <typename TIn, typename G>
    ProbabilityVector sampleDirichletDistribution(const TIn & params, G & generator) {
        ProbabilityVector retval(params.size());
 
        sampleDirichletDistribution(params, generator, retval);
 
        return retval;
    }
 
    template <typename TIn, typename TOut, typename G>
    void sampleDirichletDistribution(const TIn & params, G & generator, TOut && out) {
        assert(params.size() == out.size());
 
        double sum = 0.0;
        for (size_t i = 0; i < static_cast<size_t>(params.size()); ++i) {
            std::gamma_distribution<double> dist(params[i], 1.0);
            out[i] = dist(generator);
            sum += out[i];
        }
        out /= sum;
    }
 
    template <typename G>
    ProbabilityVector makeRandomProbability(const size_t S, G & generator) {
        ProbabilityVector b(S);
        double * bData = b.data();
        // The way this works is that we're going to generate S-1 numbers in
        // [0,1], and sort them with together with an implied 0.0 and 1.0, for
        // a total of S+1 numbers.
        //
        // The output will be represented by the differences between each pair
        // of numbers, after sorting the original vector.
        //
        // The idea is basically to take a unit vector and cut it up into
        // random parts. The size of each part is the value of an entry of the
        // output.
 
        // We must set the first element to zero even if we're later
        // overwriting it. This is to avoid bugs in case the input S is one -
        // in which case we should return a vector with a single element
        // containing 1.0.
        bData[0] = 0.0;
        for ( size_t s = 0; s < S-1; ++s )
            bData[s] = probabilityDistribution(generator);
 
        // Sort all but the implied last 1.0 which we'll add later.
        std::sort(bData, bData + S - 1);
 
        // For each number, keep track of what was in the vector there, and
        // transform it into the difference with its predecessor.
        double helper1 = bData[0], helper2;
        for ( size_t s = 1; s < S - 1; ++s ) {
            helper2 = bData[s];
            bData[s] -= helper1;
            helper1 = helper2;
        }
        // The final one is computed with respect to the overall sum of 1.0.
        bData[S-1] = 1.0 - helper1;
 
        return b;
    }
 
    inline bool checkEqualProbability(const ProbabilityVector & lhs, const ProbabilityVector & rhs) {
        const auto size = lhs.size();
        for (auto i = 0; i < size; ++i)
            if (!checkEqualSmall(lhs[i], rhs[i]))
                return false;
        return true;
    }
 
    inline double getEntropy(const ProbabilityVector & v) {
        return (v.array() * v.array().log()).sum();
    }
 
    inline double getEntropyBase2(const ProbabilityVector & v) {
        double entropy = 0.0;
        for (auto i = 0; i < v.size(); ++i)
            entropy += v[i] * std::log2(v[i]);
        return entropy;
    }
 
    ProbabilityVector projectToProbability(const Vector & v);
 
    class VoseAliasSampler {
        public:
            VoseAliasSampler(const ProbabilityVector & p);
 
            template <typename G>
            size_t sampleProbability(G & generator) const {
                const auto x = sampleDistribution_(generator);
                const int i = x;
                const auto y = x - i;
 
                if (y < prob_[i]) return i;
                return alias_[i];
            }
 
        private:
            Vector prob_;
            std::vector<size_t> alias_;
            mutable std::uniform_real_distribution<double> sampleDistribution_;
    };
}
 
#endif