sampling.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_SAMPLING_H
#define THEORETICA_SAMPLING_H
 
 
#include "../core/function.h"
#include "./prng.h"
 
 
namespace theoretica {
 
 
    using pdf_sampling_function = real(*)(const std::vector<real>&, PRNG&);
 
 
    inline real rand_uniform(real a, real b, PRNG& g, uint64_t prec = STATISTICS_RAND_PREC) {
 
        // Generate a uniform random real number in [0, 1]
        real x = (g() % prec) / static_cast<real>(prec);
 
        // Transform to target interval
        return a + (b - a) * x;
    }
 
 
    inline real rand_uniform(const std::vector<real>& theta, PRNG& g) {
 
        if(theta.size() != 2) {
            TH_MATH_ERROR("rand_uniform", theta.size(), INVALID_ARGUMENT);
            return nan();
        }
 
        return rand_uniform(theta[0], theta[1], g);
    }
 
 
    inline real rand_cointoss(PRNG& g) {
 
        return (g() % 2 == 0) ? 1 : -1;
    }
 
 
    inline real rand_cointoss(const std::vector<real>& theta, PRNG& g) {
        return rand_cointoss(g);
    }
 
 
    inline real rand_diceroll(unsigned int faces, PRNG& g) {
 
        if(faces == 0) {
            TH_MATH_ERROR("rand_diceroll", faces, INVALID_ARGUMENT);
            return nan();
        }
 
        return (g() % faces) + 1;
    }
 
 
    inline real rand_diceroll(const std::vector<real>& theta, PRNG& g) {
 
        if(theta.size() != 1) {
            TH_MATH_ERROR("rand_diceroll", theta.size(), INVALID_ARGUMENT);
            return nan();
        }
 
        if(theta[0] < 0) {
            TH_MATH_ERROR("rand_diceroll", theta[0], INVALID_ARGUMENT);
            return nan();
        }
 
        return rand_diceroll(theta[0], g);
    }
 
 
    inline real rand_trycatch(stat_function f,
        const vec<real>& theta,
        real x1, real x2,
        real y1, real y2, PRNG& g,
        unsigned int max_iter = STATISTICS_TRYANDCATCH_ITER) {
 
        real x;
        real y;
 
        unsigned int iter = 0;
 
        do {
            x = rand_uniform(x1, x2, g);
            y = rand_uniform(y1, y2, g);
            iter++;
        } while(y > f(x, theta) && iter <= max_iter);
 
        if(iter > max_iter) {
            TH_MATH_ERROR("rand_dist_tac", iter, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    inline real rand_rejectsamp(
        stat_function f, const vec<real>& theta,
        real_function p, real_function Pinv,
        PRNG& g, unsigned int max_tries = 100) {
 
        for (unsigned int i = 0; i < max_tries; ++i) {
 
            // Generate a random number following
            // the p(x) probability distribution
            // by the inverse cumulative distribution function
            const real u_1 = rand_uniform(0, 1, g);
            const real x_p = Pinv(u_1);
 
            const real u_2 = rand_uniform(0, 1, g);
 
            // Accept the sample if f(x_p)/g(x_p) > u_2
            if(u_2 * p(x_p) < f(x_p, theta))
                return x_p;
        }
 
        TH_MATH_ERROR("rand_reject_sample", max_tries, NO_ALGO_CONVERGENCE);
        return nan();
    }
 
 
    inline real rand_gaussian_polar(real mean, real sigma, PRNG& g) {
 
        static real spare;
        static bool has_spare = false;
 
        if(has_spare) {
            has_spare = false;
            return mean + spare * sigma;
        }
 
        real x, y, s;
 
        // Generate a random point inside the unit circle
        do {
 
            x = rand_uniform(-1, 1, g);
            y = rand_uniform(-1, 1, g);
            s = square(x) + square(y);
 
        } while(s >= 1 || s <= MACH_EPSILON);
 
        // Project the point
        s = sqrt(-2 * ln(s) / s);
 
        // Keep the second generated value for future calls
        spare = y * s;
        has_spare = true;
 
        return mean + sigma * x * s;
    }
 
 
    inline real rand_gaussian_boxmuller(real mean, real sigma, PRNG& g) {
 
        static real spare;
        static bool has_spare = false;
 
        if(has_spare) {
            has_spare = false;
            return mean + spare * sigma;
        }
 
        // Generate a random point inside the unit circle
        
        const real x = rand_uniform(0, 1, g);
        const real y = rand_uniform(0, 1, g);
 
        const real x_transf = sqrt(-2 * ln(x));
 
        const real u = x_transf * cos(TAU * y);
        const real v = x_transf * sin(TAU * y);
 
        spare = v;
        has_spare = true;
 
        return mean + sigma * u;
    }
 
 
    inline real rand_gaussian_clt(real mean, real sigma, PRNG& g) {
 
        // Fixed N = 12
        constexpr unsigned int N = 12;
 
        real s = 0;
        for (unsigned int i = 0; i < N; ++i)
            s += rand_uniform(-1, 1, g);
 
        // f(u) = 1/2 (in [-1, 1])
        // E[u] = 0
        // sqrt(V[u]) = 1 / sqrt(3N) = 1 / 6
 
        return mean + (s / static_cast<real>(N)) * sigma * 6;
    }
 
 
    inline real rand_gaussian_clt(
        real mean, real sigma,
        PRNG& g, unsigned int N) {
 
        real s = 0;
        for (unsigned int i = 0; i < N; ++i)
            s += rand_uniform(-1, 1, g);
 
        // f(u) = 1/2 (in [-1, 1])
        // E[u] = 0
        // sqrt(V[u]) = 1 / sqrt(3N)
 
        return mean + (s / static_cast<real>(N)) * sigma * sqrt(3 * N);
    }
 
 
    inline real rand_gaussian(real mean, real sigma, PRNG& g) {
        return rand_gaussian_polar(mean, sigma, g);
    }
 
 
    inline real rand_gaussian(const std::vector<real>& theta, PRNG& g) {
 
        if(theta.size() != 2) {
            TH_MATH_ERROR("rand_gaussian", theta.size(), INVALID_ARGUMENT);
            return nan();
        }
 
        return rand_gaussian(theta[0], theta[1], g);
    }
 
 
    inline real rand_exponential(real lambda, PRNG& g) {
 
        if(abs(lambda) < MACH_EPSILON) {
            TH_MATH_ERROR("rand_exponential", lambda, DIV_BY_ZERO);
            return nan();
        }
 
        return -ln(1 - rand_uniform(0, 1, g)) / lambda;
    }
 
 
    inline real rand_exponential(const std::vector<real>& theta, PRNG& g) {
 
        if(theta.size() != 1) {
            TH_MATH_ERROR("rand_exponential", theta.size(), INVALID_ARGUMENT);
            return nan();
        }
 
        return rand_exponential(theta[0], g);
    }
 
 
    inline real rand_rayleigh(real sigma, PRNG& g) {
 
        return sigma * sqrt(-2 * ln(1 - rand_uniform(0, 1, g)));
    }
 
 
    inline real rand_rayleigh(const std::vector<real>& theta, PRNG& g) {
 
        if(theta.size() != 1) {
            TH_MATH_ERROR("rand_rayleigh", theta.size(), INVALID_ARGUMENT);
            return nan();
        }
 
        return rand_rayleigh(theta[0], g);
    }
 
 
    inline real rand_cauchy(real mu, real alpha, PRNG& g) {
 
        return alpha * tan(PI * (rand_uniform(0, 1, g) - 0.5)) + mu;
    }
 
 
    inline real rand_cauchy(const std::vector<real>& theta, PRNG& g) {
 
        if(theta.size() != 2) {
            TH_MATH_ERROR("rand_cauchy", theta.size(), INVALID_ARGUMENT);
            return nan();
        }
 
        return rand_cauchy(theta[0], theta[1], g);
    }
 
 
    inline real rand_laplace(real mu, real b, PRNG& g) {
 
        const real u = rand_uniform(0, 0.5, g);
        return mu - b * rand_cointoss(g) * ln(1 - 2 * abs(u));
    }
 
 
    inline real rand_laplace(const std::vector<real>& theta, PRNG& g) {
 
        if(theta.size() != 2) {
            TH_MATH_ERROR("rand_laplace", theta.size(), INVALID_ARGUMENT);
            return nan();
        }
        
        return rand_laplace(theta[0], theta[1], g);
    }
 
 
    inline real rand_pareto(real x_m, real alpha, PRNG& g) {
 
        return x_m / powf(1 - rand_uniform(0, 1, g), 1.0 / alpha);
    }
 
 
    inline real rand_pareto(const std::vector<real>& theta, PRNG& g) {
 
        if(theta.size() != 2) {
            TH_MATH_ERROR("rand_pareto", theta.size(), INVALID_ARGUMENT);
            return nan();
        }
 
        return rand_pareto(theta[0], theta[1], g);
    }
 
 
    struct pdf_sampler {
 
        pdf_sampling_function pdf;
 
        std::vector<real> theta;
 
        PRNG& generator;
 
 
        pdf_sampler(
            pdf_sampling_function pdf,
            const std::vector<real>& theta,
            PRNG& generator) : pdf(pdf), theta(theta), generator(generator) {}
 
 
        inline real next() {
            return pdf(theta, generator);
        }
 
        inline real operator()() {
            return next();
        }
 
        template<typename Vector>
        inline void fill(Vector& x, size_t N) {
 
            for (size_t i = 0; i < N; ++i)
                x[i] = next();
        }
 
 
        template<typename Vector>
        inline void fill(Vector& x) {
 
            for (size_t i = 0; i < x.size(); ++i)
                x[i] = next();
        }
 
 
        inline pdf_sampler& operator>>(real& x) {
            x = next();
            return *this;
        }
 
 
        static pdf_sampler uniform(real a, real b, PRNG& generator) {
            return pdf_sampler(rand_uniform, {a, b}, generator);
        }
 
 
        static pdf_sampler gaussian(real mean, real sigma, PRNG& generator) {
            return pdf_sampler(rand_gaussian, {mean, sigma}, generator);
        }
 
 
        static pdf_sampler exponential(real lambda, PRNG& generator) {
            return pdf_sampler(rand_exponential, {lambda}, generator);
        }
 
 
        static pdf_sampler cauchy(real mu, real alpha, PRNG& generator) {
            return pdf_sampler(rand_cauchy, {mu, alpha}, generator);
        }
 
 
        static pdf_sampler rayleigh(real sigma, PRNG& generator) {
            return pdf_sampler(rand_rayleigh, {sigma}, generator);
        }
 
 
        static pdf_sampler pareto(real x_m, real alpha, PRNG& generator) {
            return pdf_sampler(rand_pareto, {x_m, alpha}, generator);
        }
 
 
        static pdf_sampler laplace(real mu, real b, PRNG& generator) {
            return pdf_sampler(rand_laplace, {mu, b}, generator);
        }
 
    };
 
    inline real metropolis(
        real_function pdf, pdf_sampler& g,
        real x0, PRNG& rnd, unsigned int depth = STATISTICS_METROPOLIS_DEPTH) {
 
        real current = x0, next;
 
        for(unsigned int i = 0; i < depth; i++) {
            
            // Computes the next step
            next = current + g();
 
            // Checks acceptance rate
            if(rand_uniform(0, 1, rnd) * pdf(current) <= pdf(next))
                current = next;
        }
 
        return current;
    }
 
 
    inline real metropolis(real_function pdf, pdf_sampler& g,
        real x0, unsigned int depth = STATISTICS_METROPOLIS_DEPTH) {
        return metropolis(pdf, g, x0, g.generator, depth);
    }
 
}
 
#endif