prob_distribution.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_distribution.h 18558 2015-02-13 21:13:12Z jguan1 $ */

#ifndef PROB_DISTRIBUTION_H_INCLUDED
#define PROB_DISTRIBUTION_H_INCLUDED

#include <boost/version.hpp>
#include <boost/math/distributions.hpp>
#include <boost/math/distributions/inverse_gamma.hpp>
#include <boost/math/distributions/negative_binomial.hpp>
#include <map>
#include <cmath>
#include <vector>
#include "prob_cpp/prob_util.h"
#include "m_cpp/m_vector.h"
#include "m_cpp/m_vector_d.h"
#include "m_cpp/m_matrix.h"
#include "n_cpp/n_cholesky.h"
#include "l_cpp/l_exception.h"
#include "l_cpp/l_std_parallel.h"

#include "g/g_area.h"
#include "g_cpp/g_quaternion.h"

namespace kjb{

/*============================================================================
  These distributions are part of the boost-math library, defined in
  boost/math/distributions.hpp. See boost's documentation for details.
  They are typdef'd to match the naming scheme of the KJB
  ----------------------------------------------------------------------------*/

typedef boost::math::bernoulli              Bernoulli_distribution;
typedef boost::math::beta_distribution<>    Beta_distribution;
typedef boost::math::binomial               Binomial_distribution;
typedef boost::math::chi_squared            Chi_square_distribution;
typedef boost::math::exponential            Exponential_distribution;
typedef boost::math::gamma_distribution<>   Gamma_distribution;
typedef boost::math::normal                 Gaussian_distribution;
typedef boost::math::laplace                Laplace_distribution;
typedef boost::math::normal                 Normal_distribution;
typedef boost::math::poisson                Poisson_distribution;
typedef boost::math::uniform                Uniform_distribution;
typedef boost::math::inverse_gamma_distribution<>   Inverse_gamma_distribution;

// standard normal distribution
extern const Gaussian_distribution STD_NORMAL;

#if BOOST_VERSION >= 104600
typedef boost::math::geometric_distribution<double> Geometric_distribution;
#else
class Geometric_distribution : public boost::math::negative_binomial
{
public:
    Geometric_distribution(double p) : 
        boost::math::negative_binomial(1, p)
    { }
};
#endif

/* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
class Beta_binomial_distribution
{
public:
    Beta_binomial_distribution(size_t n_, double alpha_, double beta_) : 
        n(n_), alpha(alpha_), beta(beta_)
    { }

    double n, alpha, beta;
};

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

template<class Distribution>
struct Distribution_traits
{
    typedef double type;
};

/*============================================================================
  The following are univariate distributions that are not present in the
  boost library.
  ----------------------------------------------------------------------------*/

template<class T = int>
class Categorical_distribution
{
    typedef Categorical_distribution Self;
    typedef std::pair<const T*, double> Cdf_pair;
    typedef typename std::map<T, double>::value_type Map_pair;

public:
    Categorical_distribution() :
        db(),
        cdf_(),
        total_weight_(0)
    {}

    Categorical_distribution(const std::vector<double>& ps, const T& t1 = T(1));

    Categorical_distribution(const Vector& ps, const T& = T(1));

    Categorical_distribution
    (
        const std::vector<T>& values,
        const std::vector<double>& probabilities
    );

    /*template <class Value_iterator, class Prob_iterator>
    Categorical_distribution
    (
        Value_iterator values_begin,
        Value_iterator values_end,
        Prob_iterator probabilities_begin
    );*/

    Categorical_distribution(const std::map<T, double>& d) : db(d)
    {
        init_cdf_();
    }

    Categorical_distribution(const std::map<T, size_t>& m)
    {
        for(typename std::map<T,size_t>::const_iterator it = m.begin();
            it != m.end(); ++it)
        {
            db[it->first] = static_cast<double>(it->second);
        }
        init_cdf_();
    }
    
    Categorical_distribution(const T& min, const T& max, const T& step);

    Categorical_distribution(const T& x1, const T& x2, double p1, double p2)
    {
        db[x1] = p1;
        db[x2] = p2;
        init_cdf_();
    }

    Categorical_distribution(const Categorical_distribution<T>& cd) :
        db(cd.db)
    {
        init_cdf_();
    }

    Categorical_distribution<T>& operator=(const Categorical_distribution<T>& cd)
    {
        if(this != &cd)
        {
            db = cd.db;
            init_cdf_();
        }

        return *this;
    }

    static Categorical_distribution<T> construct_from_log_map(
        const std::map<T, double>& log_map
        )
    {
        typedef typename std::map<T, double>::value_type Pair;
        Categorical_distribution<T> result;

        result.db = log_map;

        // get the sum of the log-weights, without leaving log space, which mitigates precision loss
        // when weights are very small.

        std::vector<double> log_probs;
        log_probs.reserve(log_map.size());

        typedef typename std::map<T,double>::const_iterator iterator;
        for(iterator it = log_map.begin(); it != log_map.end(); ++it)
        {
            log_probs.push_back(it->second);
        }

        double log_total = kjb::log_sum(log_probs.begin(), log_probs.end());

        // normalize by subtracting log-total.  then exponentiate to get out of log-space
        for(iterator it = result.db.begin(); it != result.db.end(); ++it)
            it->second = exp(it->second - log_total);

        // TODO: init_cdf_ re-normalizes here; would be nice to factor the normalization
        // bits out of init_cdf_.  Note: we can't avoid normalizing above, because if we don't
        // we risk losing precision when exponentiating.
        result.init_cdf_();
        return result;
    }

    size_t size() const
    {
        return db.size();
    }

    void erase(const T& key)
    {
        size_t count = db.erase(key);
        if(count == 0) 
            KJB_THROW_2(Illegal_argument, "Key not found");
        init_cdf_();
    }

    void remove(const T& key)
    {
        erase(key);
    }

    void insert (const T& key, double weight = 1.0)
    {
        size_t count = db.count(key);
        if(count > 0)
            KJB_THROW_2(Illegal_argument, "Key collision");

        const T* stored_key = &db.insert(std::make_pair(key, weight)).first->first;

        if(cdf_.empty())
            total_weight_ = weight;
        else
            total_weight_ += weight;

        cdf_.push_back(Cdf_pair(stored_key, total_weight_));
    }

    void add(const T& key, double weight = 1.0)
    {
        insert(key, weight);
    }

    ~Categorical_distribution()
    {}

    // friends functions; these need to be friends because they
    // access the private parts
    template<class U>
    friend double pdf(const Categorical_distribution<U>& dist, const U& x);

    template<class U>
    friend double cdf(const Categorical_distribution<U>& dist, const U& x);

    template<class U>
    friend U sample(const Categorical_distribution<U>& dist);

private:

    void init_cdf_()
    {
        using namespace boost;

        cdf_.resize(db.size());
        if(cdf_.size() == 0) return;

        //typedef typename std::map<T, double>::const_iterator Const_iterator;

        // copy into cdf
        kjb_parallel_std::transform(
            db.begin(),
            db.end(),
            cdf_.begin(),
            make_cdf_pair()
        );

        // convert pdf values into cdf values
        kjb_parallel_std::partial_sum(
            cdf_.begin(),
            cdf_.end(),
            cdf_.begin(),
            partial_sum_second<Cdf_pair>()
        );

        total_weight_ = cdf_.back().second;

        // ensure that at least one element had significant mass,
        // so we can normalize the masses
        if(total_weight_ < FLT_EPSILON)
        {
            KJB_THROW_2(Illegal_argument, "Items have zero probability mass.");
        }


        // normalize to 1.0
        // update: we now store CDF unnormalized, to make dynamic
        // adding of elements easier.
//        kjb_parallel_std::for_each(
//            cdf_.begin(),
//            cdf_.end(),
//            multiply_second_by_scalar<Cdf_pair>(1.0/total_weight_)
//        );

        // update: we now store elements unnormalized, to make dynamic
        // adding of elements easier.
//        // normalize pdf map to 1.0
//        kjb_parallel_std::for_each(
//            db.begin(),
//            db.end(),
//            multiply_second_by_scalar<Map_pair>(1.0/cdf_.back().second)
//        );
    }


    // utility functors
    struct make_cdf_pair : public std::unary_function<Cdf_pair, const Map_pair&>
    {
        Cdf_pair operator()(const Map_pair& pair) const
        {
            return std::make_pair(&pair.first, pair.second);
        }
    };

    struct cdf_pair_less_than_scalar
    {
        bool operator()(const Cdf_pair& op1, double op2) const
        {
            return op1.second < op2;
        }
    };

    struct cdf_key_less_than_scalar
    {
        bool operator()(const Cdf_pair& op1, const T& op2) const
        {
            return op1.first < op2;
        }
    };

    template <class U>
    struct multiply_second_by_scalar
    {
        multiply_second_by_scalar(double scalar) :
            scalar_(scalar) 
        { }

        void operator()(U& in) const
        {
            in.second *= scalar_;
        }
        double scalar_;
    };


    template <class U>
    struct partial_sum_second : public std::binary_function<U, const U&, const U&>
    {
        U operator()(const U& op1, const U& op2)
        {
            U result = op2;
            result.second += op1.second;
            return result;
        }
    };
private:
    std::map<T, double> db;
    std::vector<Cdf_pair> cdf_;
    double total_weight_;
};

/* non-inline member functions */

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

template<class T>
Categorical_distribution<T>::Categorical_distribution
(
    const std::vector<double>& ps,
    const T& t1
)
{
    for(std::size_t i = 0; i < ps.size(); i++)
    {
        db[static_cast<T>(i + t1)] = ps[i];
    }
    init_cdf_();
}

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

template<class T>
Categorical_distribution<T>::Categorical_distribution
(
    const Vector& ps,
    const T& t1
)
{
    for(int i = 0; i < ps.get_length(); i++)
    {
        db[static_cast<T>(i + t1)] = ps[i];
    }
    init_cdf_();
}


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

template <class T>
Categorical_distribution<T>::Categorical_distribution
(
    const std::vector<T>& values,
    const std::vector<double>& probabilities)
{
    for(std::size_t i = 0; i < probabilities.size(); i++)
    {
        db[values[i]] = probabilities[i];
    }
    init_cdf_();
}

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

/*template <class T>
template <class Value_iterator, class Prob_iterator>
Categorical_distribution<T>::Categorical_distribution
(
    Value_iterator values_begin,
    Value_iterator values_end,
    Prob_iterator probabilities_begin)
{
    Prob_iterator prob_it = probabilities_begin;
    for(Value_iterator it = values_begin; it != values_end; ++it, ++prob_it)
    {
        db[*it] = *prob_it; 
    }
    init_cdf_();
}*/

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

template<class T>
Categorical_distribution<T>::Categorical_distribution
(
    const T& min,
    const T& max,
    const T& step
)
{
    double p = 1.0 / static_cast<int>((max - min + 1) / step);
    for(T x = min; x <= max; x += step)
    {
        db[x] = p;
    }
    init_cdf_();
}


template<class T>
struct Distribution_traits<Categorical_distribution<T> >
{
    typedef T type;
};

/*============================================================================
  The following are multivariate distributions that are not present in the
  boost library.
  ----------------------------------------------------------------------------*/

// forward declarations -- ignore
class MV_gaussian_distribution;
template<class T, class G, class F> class Conditional_distribution;
class MV_normal_on_normal_dependence;
typedef Conditional_distribution<MV_gaussian_distribution, MV_gaussian_distribution, MV_normal_on_normal_dependence> MV_gaussian_conditional_distribution;

class MV_gaussian_distribution
{
public:

    MV_gaussian_distribution(int d) : 
        mean(d, 0.0),
        cov_mat(create_identity_matrix(d)),
        type(INDEPENDENT),
        abs_det(1.0),
        log_abs_det(0.0)
    {}

    MV_gaussian_distribution(const Vector& mu, const Matrix& sigma) :
        mean(mu),
        cov_mat(sigma),
        type(NORMAL),
        abs_det(-1.0),
        log_abs_det(-std::numeric_limits<double>::max())
    {
        IFT(mean.get_length() == cov_mat.get_num_rows()
                && cov_mat.get_num_rows() == cov_mat.get_num_cols(),
            Dimension_mismatch,
            "Mean vector dimension must equal covariance matrix's"
                "number of rows and columns.");
    }

    MV_gaussian_distribution(const Vector& mu, const Vector& sigma_diag) :
        mean(mu),
        cov_mat(create_diagonal_matrix(sigma_diag)),
        type(INDEPENDENT),
        abs_det(-1.0),
        log_abs_det(-std::numeric_limits<double>::max())
    {
        IFT(mean.get_length() == cov_mat.get_num_rows(), Dimension_mismatch,
            "Mean vector dimension must equal covariance matrix's"
                "number of rows and columns.");
    }

    const Vector& get_mean() const 
    {
        return mean;
    }

    const Matrix& get_covariance_matrix() const 
    {
        return cov_mat;
    }

    void set_mean(const Vector& m)
    {
        IFT(cov_mat.get_num_rows() == m.get_length(), Illegal_argument,
            "mean vector must have same size as covariance matrix.");

        mean = m;
    }

    void set_covariance_matrix(const Matrix& S)
    {
        IFT(S.get_num_rows() == mean.get_length(), Illegal_argument,
            "mean vector must have same size as covariance matrix.");

        type = NORMAL;

        cov_mat = S;
        cov_inv.resize(0, 0);
        cov_chol.resize(0, 0);
        abs_det = -1;
        log_abs_det = -std::numeric_limits<double>::max();
    }

    void set_covariance_matrix(const Matrix& S, const Vector& m)
    {
        IFT(S.get_num_rows() == m.get_length(), Illegal_argument,
            "mean vector must have same size as covariance matrix.");

        type = NORMAL;

        mean = m;
        cov_mat = S;
        cov_inv.resize(0, 0);
        cov_chol.resize(0, 0);
        abs_det = -1;
        log_abs_det = -std::numeric_limits<double>::max();
    }

    int get_dimension() const
    {
        return mean.get_length();
    }

    MV_gaussian_conditional_distribution conditional(int i) const;

private:
    void update_cov_inv() const
    {
        if(cov_inv.get_num_rows() == 0)
        {
            cov_inv = matrix_inverse(cov_mat);
        }
    }

    void update_cov_chol() const
    {
        if(cov_chol.get_num_rows() == 0)
        {
            cov_chol = cholesky_decomposition(cov_mat);
        }
    }

    void update_abs_det() const
    {
        if(abs_det == -1.0)
        {
            abs_det = cov_mat.abs_of_determinant();
        }
    }

    void update_log_abs_det() const;

    friend double pdf(const MV_gaussian_distribution&, const Vector&);

    friend double log_pdf(const MV_gaussian_distribution&, const Vector&);

    friend Vector sample(const MV_gaussian_distribution&);

private:
    Vector mean;
    Matrix cov_mat;
    int type;
    mutable Matrix cov_inv;
    mutable Matrix cov_chol;
    mutable double abs_det;
    mutable double log_abs_det;

    const static int NORMAL = 0;
    const static int INDEPENDENT = 1;
};

//We should also be able to call it normal distribution
typedef MV_gaussian_distribution MV_normal_distribution;

template<>
struct Distribution_traits<MV_gaussian_distribution>
{
    typedef Vector type;
};

class Log_normal_distribution
{
public:
    Log_normal_distribution() : 
        dist_(0, 1)
    {}

    Log_normal_distribution(double mu, double sigma) :
        dist_(mu, sigma)
    {}

    friend double pdf(const Log_normal_distribution&, double);

    friend double cdf(const Log_normal_distribution&, double);

    friend double sample(const Log_normal_distribution&);
private:
    Normal_distribution dist_;
};

template<>
struct Distribution_traits<Log_normal_distribution>
{
    typedef double type;
};


/*============================================================================
  The following are other distributions (e.g., mixture).
  ============================================================================*/

template<class Distribution>
class Mixture_distribution
{
private:
    // a convenient typedef
    typedef typename Distribution_traits<Distribution>::type type;

public:

    Mixture_distribution
    (
        const std::vector<Distribution>& distributions,
        const std::vector<double>& coefficients
    ) :
        dists(distributions), coeffs(coefficients)
    {}

    Mixture_distribution(const std::vector<Distribution>& distributions) :
        dists(distributions), coeffs(1.0, distributions.size(), 1.0)
    {}

    Mixture_distribution
    (
        const Distribution& dist1,
        const Distribution& dist2,
        double pi1
    ) :
        coeffs(1, 2, pi1, 1 - pi1)
    {
        dists.push_back(dist1);
        dists.push_back(dist2);
    }

    // friends functions; these need to be friends because they access
    // the private parts (lol)
    template<class Dist>
    friend
    double pdf
    (
        const kjb::Mixture_distribution<Dist>& dist,
        const typename Distribution_traits<Mixture_distribution<Dist> >::type& x
    );

    template<class Dist>
    friend
    double cdf
    (
        const kjb::Mixture_distribution<Dist>& dist,
        const typename Distribution_traits<Mixture_distribution<Dist> >::type& x
    );

    template<class Dist>
    friend
    typename Distribution_traits<Mixture_distribution<Dist> >::type
        sample(const kjb::Mixture_distribution<Dist>& dist);

private:
    std::vector<Distribution> dists;
    Categorical_distribution<size_t> coeffs;
};

template<class Distribution>
struct Distribution_traits<Mixture_distribution<Distribution> >
{
    typedef typename Distribution_traits<Distribution>::type type;
};

/* /////////////////////////////////////////////////////////////////////////// */

template <size_t D>
class Uniform_sphere_distribution
{
public:
    Uniform_sphere_distribution()  {}
};

template<size_t D>
struct Distribution_traits<Uniform_sphere_distribution<D> >
{
    typedef typename kjb::Vector_d<D> type;
};

/* /////////////////////////////////////////////////////////////////////////// */

template <size_t D>
class Von_mises_fisher_distribution
{
public:
    Von_mises_fisher_distribution(const kjb::Vector_d<D>& mean, double kappa) :
        mu_(mean.normalized()),
        kappa_(kappa)
    {}

    const kjb::Vector_d<D>& mu() const { return mu_; }

    double kappa() const { return kappa_; }

private:
    kjb::Vector_d<D> mu_;
    double kappa_;
}; // Von_mises_fisher_distribution

template<size_t D>
struct Distribution_traits<Von_mises_fisher_distribution<D> >
{
    typedef typename kjb::Vector_d<D> type;
};

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

class Chinese_restaurant_process
{
public:
    typedef std::vector<std::vector<size_t> > Type;

public:
    Chinese_restaurant_process(double concentration, size_t nc) :
        theta_(concentration), nc_(nc)
    {}

    double concentration() const { return theta_; }
    size_t num_customers() const { return nc_; }

private:
    double theta_;
    double nc_;
};

template<>
struct Distribution_traits<Chinese_restaurant_process>
{
    typedef Chinese_restaurant_process::Type type;
};

// useful typedef
typedef Chinese_restaurant_process Crp;

} //namespace kjb


#endif /*PROB_DISTRIBUTION_H_INCLUDED */