prob_sample.h

Go to the documentation of this file.

/* ======================================================================== *
 |                                                                          |
 | Copyright (c) 2007-2010, by members of University of Arizona Computer    |
 | Vision group (the authors) including                                     |
 |       Kobus Barnard, Kyle Simek, Ernesto Brau.                           |
 |                                                                          |
 | For use outside the University of Arizona Computer Vision group please   |
 | contact Kobus Barnard.                                                   |
 |                                                                          |
 * ======================================================================== */

/* $Id: prob_sample.h 18761 2015-04-04 00:47:45Z cdawson $ */

#ifndef PROB_SAMPLE_H_INCLUDED
#define PROB_SAMPLE_H_INCLUDED

#include <prob_cpp/prob_distribution.h>
#include <prob_cpp/prob_pdf.h>
#include <m_cpp/m_vector.h>
#include <m_cpp/m_vector_d.h>
#include <m_cpp/m_matrix_d.h>
#include <g_cpp/g_util.h>
#include <l_cpp/l_util.h>
#include <boost/random.hpp>
#include <algorithm>
#include <vector>

namespace kjb {

extern const unsigned int DEFAULT_SEED;

// Basic random generator -- used for sampling -- belongs to the cpp file
//typedef boost::minstd_rand Base_generator_type;
typedef boost::mt19937 Base_generator_type;
extern Base_generator_type basic_rnd_gen;
//extern boost::uniform_01<Base_generator_type> uni01;

/* /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ */

inline void seed_sampling_rand(unsigned int x0)
{
    basic_rnd_gen.seed(x0);
//    uni01 = boost::uniform_01<Base_generator_type>(basic_rnd_gen);
}

inline double sample(const Uniform_distribution& dist)
{
    typedef boost::uniform_real<> Distribution_type;
    typedef boost::variate_generator<Base_generator_type&, Distribution_type> Rng;
    Rng rng(basic_rnd_gen, Distribution_type(dist.lower(), dist.upper()));
    return rng();
}

template<class Distro>
inline double sample(const Distro& dist)
{
    // this sometimes results in infinite recursion inside of boost.  No ideas :-(
    double u = sample(Uniform_distribution());
    return quantile(dist, u);
}

inline double sample(const Bernoulli_distribution& dist)
{
    typedef boost::bernoulli_distribution<> Distribution_type;
    typedef boost::variate_generator<Base_generator_type&, Distribution_type> Rng;

    Rng rng(basic_rnd_gen, Distribution_type(dist.success_fraction()));
    return rng();
}

inline double sample(const Binomial_distribution& dist)
{
    typedef boost::binomial_distribution<> Distribution_type;
    typedef boost::variate_generator<Base_generator_type&, Distribution_type> Rng;

    Rng rng(basic_rnd_gen, Distribution_type(dist.trials(), dist.success_fraction()));
    return rng();
}

inline double sample(const Exponential_distribution& dist)
{
    typedef boost::exponential_distribution<> Distribution_type;
    typedef boost::variate_generator<Base_generator_type&, Distribution_type> Rng;

    Rng rng(basic_rnd_gen, Distribution_type(dist.lambda()));
    return rng();
}

inline double sample(const Gaussian_distribution& dist)
{
    typedef boost::normal_distribution<> Distribution_type;
    typedef boost::variate_generator<Base_generator_type&, Distribution_type> Rng;

    Rng rng(basic_rnd_gen, Distribution_type(dist.mean(), dist.standard_deviation()));
    return rng();
}

inline double sample(const Poisson_distribution& dist)
{
    typedef boost::poisson_distribution<> Distribution_type;
    typedef boost::variate_generator<Base_generator_type&, Distribution_type> Rng;

    Rng rng(basic_rnd_gen, Distribution_type(dist.mean()));
    return rng();
}

#if BOOST_VERSION < 103400
inline double sample(const Geometric_distribution& dist)
{
#warning "sample(Geometric_distribution) not implemented in Boost prior to 1.034"
    KJB_THROW(Not_implemented);
}
#else
inline double sample(const Geometric_distribution& dist)
{
    // this might not be available in older boost versions but exactly which version it was added is kind of a pain.  
    // If you can't compile because this isn't available, look in /usr/local/boost/version.hpp and determine which version you have.
#if BOOST_VERSION >= 104700
    double dist_param = dist.success_fraction();
#else /*  BOOST_VERSION >= 103400 */
    KJB(UNTESTED_CODE());
    double dist_param = 1.0-dist.success_fraction();
#endif

    typedef boost::geometric_distribution<> Distribution_type;
    typedef boost::variate_generator<Base_generator_type&, Distribution_type> Rng;

    Rng rng(basic_rnd_gen, Distribution_type(dist_param));
    return rng();
}
#endif /* BOOST_VERSION < 10340 */

template<class T>
T sample(const Categorical_distribution<T>& dist)
{
    // TODO: This is currently O(log n), but approaches exist that
    // make this constant, see the Alias Table methods described here:
    // http://www.keithschwarz.com/darts-dice-coins/
    // Also boost's discrete_distribution implements most of what
    // this class does, so we could replace this with that.
    double u = sample(Uniform_distribution());

    using namespace boost;

    if(dist.cdf_.empty())
    {
        KJB_THROW_2(Runtime_error, "Attempting to sample from an empty set.");
    }

    assert(dist.cdf_.size() == dist.db.size()); // init_cdf() was called
    double total_weight = dist.cdf_.back().second;
    u *= total_weight;

    typedef Categorical_distribution<T> Cd;
    typedef typename Cd::Cdf_pair Cdf_pair;
    typedef typename std::vector<Cdf_pair>::const_iterator Const_iterator;

    // use binary search to find first element in CDF that exceeds u
    Const_iterator r = 
        std::lower_bound(
            dist.cdf_.begin(),
            dist.cdf_.end(),
            u,
            typename Cd::cdf_pair_less_than_scalar());

    assert(r != dist.cdf_.end());

    return *r->first;

}

template<class Iterator, class distance_type>
inline
Iterator element_uar(Iterator first, distance_type size)
{
    int dist = kjb::sample(kjb::Categorical_distribution<size_t>(0, size - 1, 1));
    Iterator p = first;
    std::advance(p, dist);

    return p;
}

Vector sample(const MV_gaussian_distribution& dist);

inline double sample(const Log_normal_distribution& dist)
{
    return exp(sample(dist.dist_));
}

template<class Distribution>
inline
typename Distribution_traits<Mixture_distribution<Distribution> >::type
    sample(const Mixture_distribution<Distribution>& dist)
{
    size_t i = sample(dist.coeffs);
    return sample(dist.dists[i - 1]);
}

template <size_t D> inline kjb::Vector_d<D> sample(const Uniform_sphere_distribution<D>& /* dist */)
{
    kjb::Vector_d<D> result;

    static const Gaussian_distribution norm_dist(0, 1);

    for(size_t d = 0; d < D; ++d)
    {
        result[d] = kjb::sample(norm_dist);
    }

    return result.normalize();
}


template <size_t D, class Vector_d_iterator>
inline 
void sample(const Von_mises_fisher_distribution<D>& dist, size_t n, Vector_d_iterator it)
{
    size_t m = D;
    static const kjb::Vector_d<D> z_hat = kjb::create_unit_vector<D>(D-1);

    const Vector_d<D>& mu_ = dist.mu();
    double kappa_ = dist.kappa();

    // TODO: Quaternion may be slower than is needed here.  Consider
    // constructing the rotation matrix by hand or move q into a member variable
    // I replaced it with a function I wrote to find rotation matrices
    // between two directions in any dimension. I haven't tested it as
    // thoroughly as it needs, but I'm pretty confident it's correct.
    //                                       -- Ernesto [2014/10/28]
    //kjb::Quaternion q; q.set_from_directions(z_hat, mu_);
    Matrix_d<D, D> R = geometry::get_rotation_matrix(z_hat, mu_);

    double b = (-2*kappa_ + sqrt(4*kappa_*kappa_ + (m-1)*(m-1)))/(m-1);
    double x0 = (1-b)/(1+b);
    double c = kappa_*x0 + (m-1)*log(1-x0*x0);

    Uniform_sphere_distribution<D-1> subdir_dist;
    for(size_t i = 0; i < n; ++i)
    {
        double W, U;
        do
        {
            // rejection sampling for height on unit sphere,
            // distributed as f(w) = k * exp(kappa*w) ; for w \in [-1, 1], zero elsewhere
            // see Sungkyu Jung's treatment:.
            // http://www.stat.pitt.edu/sungkyu/software/randvonMisesFisher3.pdf
            double Z = kjb::sample(kjb::Beta_distribution((m-1)/2,(m-1)/2));
            W = (1-(1+b)*Z)/(1-(1-b)*Z);
            U = kjb::sample(kjb::Uniform_distribution());
        } while(kappa_*W + (m-1)*log(1-x0*W) - c < log(U));


        kjb::Vector_d<D-1> subdir = kjb::sample(subdir_dist);
        subdir *= sqrt(1- W*W);

        kjb::Vector_d<D> dir;
        std::copy(subdir.begin(), subdir.end(), dir.begin());
        dir.back() = W;

        //dir = q.rotate(dir);
        dir = R*dir;
        *it = dir;
        ++it;
    }
}

template <size_t D>
Vector_d<D> sample(const Von_mises_fisher_distribution<D>& dist)
{
    Vector_d<D> result;
    sample(dist, 1, &result);
    return result;
}

Crp::Type sample(const Chinese_restaurant_process& crp);

size_t sample_occupied_tables(const Chinese_restaurant_process& crp);
    
} //namespace kjb

#endif /*PROB_SAMPLE_H_INCLUDED */