Commit bef459e7 authored by Manuel Günther's avatar Manuel Günther
Browse files

#WARNING: implementation of bob.core.lognormal was changed to new version;...

#WARNING: implementation of bob.core.lognormal was changed to new version; fixed implementations of random number distributions (added boost 1.56 headers in namespace bob::core::random); added C++ API for random number distributions in bob::core::random
parent c23a9441
include LICENSE README.rst bootstrap.py buildout.cfg requirements.txt
recursive-include doc conf.py *.rst
recursive-include bob *.cpp *.h
recursive-include bob *.cpp *.h *.hpp
......@@ -9,7 +9,7 @@
#define BOB_CORE_ARRAY_RANDOM_H
#include <blitz/array.h>
#include <boost/random.hpp>
#include <bob.core/random.h>
namespace bob { namespace core { namespace array {
/**
......@@ -29,12 +29,10 @@ template <typename T, int N>
void randn(boost::mt19937& rng, blitz::Array<T,N>& dest,
const T mean=0., const T var=1.)
{
boost::normal_distribution<T> distrib(mean, var);
boost::variate_generator<boost::mt19937, boost::normal_distribution<T> >
vg(rng, distrib);
bob::core::random::normal_distribution<T> distrib(mean, var);
typedef typename blitz::Array<T,N>::iterator BA_iterator;
for (BA_iterator it=dest.begin(); it !=dest.end(); ++it)
*it = vg();
*it = distrib(rng);
}
/**
......
/* boost random/binomial_distribution.hpp header file
*
* Copyright Steven Watanabe 2010
* Distributed under the Boost Software License, Version 1.0. (See
* accompanying file LICENSE_1_0.txt or copy at
* http://www.boost.org/LICENSE_1_0.txt)
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*/
#ifndef BOB_BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
#define BOB_BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
#include <boost/config/no_tr1/cmath.hpp>
#include <cstdlib>
#include <iosfwd>
#include <boost/random/detail/config.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/detail/disable_warnings.hpp>
namespace bob {
namespace core {
namespace random {
namespace detail {
template<class RealType>
struct binomial_table {
static const RealType table[10];
};
template<class RealType>
const RealType binomial_table<RealType>::table[10] = {
0.08106146679532726,
0.04134069595540929,
0.02767792568499834,
0.02079067210376509,
0.01664469118982119,
0.01387612882307075,
0.01189670994589177,
0.01041126526197209,
0.009255462182712733,
0.008330563433362871
};
}
/**
* The binomial distribution is an integer valued distribution with
* two parameters, @c t and @c p. The values of the distribution
* are within the range [0,t].
*
* The distribution function is
* \f$\displaystyle P(k) = {t \choose k}p^k(1-p)^{t-k}\f$.
*
* The algorithm used is the BTRD algorithm described in
*
* @blockquote
* "The generation of binomial random variates", Wolfgang Hormann,
* Journal of Statistical Computation and Simulation, Volume 46,
* Issue 1 & 2 April 1993 , pages 101 - 110
* @endblockquote
*/
template<class IntType = int, class RealType = double>
class binomial_distribution {
public:
typedef IntType result_type;
typedef RealType input_type;
class param_type {
public:
typedef binomial_distribution distribution_type;
/**
* Construct a param_type object. @c t and @c p
* are the parameters of the distribution.
*
* Requires: t >=0 && 0 <= p <= 1
*/
explicit param_type(IntType t_arg = 1, RealType p_arg = RealType (0.5))
: _t(t_arg), _p(p_arg)
{}
/** Returns the @c t parameter of the distribution. */
IntType t() const { return _t; }
/** Returns the @c p parameter of the distribution. */
RealType p() const { return _p; }
#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
/** Writes the parameters of the distribution to a @c std::ostream. */
template<class CharT, class Traits>
friend std::basic_ostream<CharT,Traits>&
operator<<(std::basic_ostream<CharT,Traits>& os,
const param_type& parm)
{
os << parm._p << " " << parm._t;
return os;
}
/** Reads the parameters of the distribution from a @c std::istream. */
template<class CharT, class Traits>
friend std::basic_istream<CharT,Traits>&
operator>>(std::basic_istream<CharT,Traits>& is, param_type& parm)
{
is >> parm._p >> std::ws >> parm._t;
return is;
}
#endif
/** Returns true if the parameters have the same values. */
friend bool operator==(const param_type& lhs, const param_type& rhs)
{
return lhs._t == rhs._t && lhs._p == rhs._p;
}
/** Returns true if the parameters have different values. */
friend bool operator!=(const param_type& lhs, const param_type& rhs)
{
return !(lhs == rhs);
}
private:
IntType _t;
RealType _p;
};
/**
* Construct a @c binomial_distribution object. @c t and @c p
* are the parameters of the distribution.
*
* Requires: t >=0 && 0 <= p <= 1
*/
explicit binomial_distribution(IntType t_arg = 1,
RealType p_arg = RealType(0.5))
: _t(t_arg), _p(p_arg)
{
init();
}
/**
* Construct an @c binomial_distribution object from the
* parameters.
*/
explicit binomial_distribution(const param_type& parm)
: _t(parm.t()), _p(parm.p())
{
init();
}
/**
* Returns a random variate distributed according to the
* binomial distribution.
*/
template<class URNG>
IntType operator()(URNG& urng) const
{
if(use_inversion()) {
if(0.5 < _p) {
return _t - invert(_t, 1-_p, urng);
} else {
return invert(_t, _p, urng);
}
} else if(0.5 < _p) {
return _t - generate(urng);
} else {
return generate(urng);
}
}
/**
* Returns a random variate distributed according to the
* binomial distribution with parameters specified by @c param.
*/
template<class URNG>
IntType operator()(URNG& urng, const param_type& parm) const
{
return binomial_distribution(parm)(urng);
}
/** Returns the @c t parameter of the distribution. */
IntType t() const { return _t; }
/** Returns the @c p parameter of the distribution. */
RealType p() const { return _p; }
/** Returns the smallest value that the distribution can produce. */
IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
/** Returns the largest value that the distribution can produce. */
IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const { return _t; }
/** Returns the parameters of the distribution. */
param_type param() const { return param_type(_t, _p); }
/** Sets parameters of the distribution. */
void param(const param_type& parm)
{
_t = parm.t();
_p = parm.p();
init();
}
/**
* Effects: Subsequent uses of the distribution do not depend
* on values produced by any engine prior to invoking reset.
*/
void reset() { }
#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
/** Writes the parameters of the distribution to a @c std::ostream. */
template<class CharT, class Traits>
friend std::basic_ostream<CharT,Traits>&
operator<<(std::basic_ostream<CharT,Traits>& os,
const binomial_distribution& bd)
{
os << bd.param();
return os;
}
/** Reads the parameters of the distribution from a @c std::istream. */
template<class CharT, class Traits>
friend std::basic_istream<CharT,Traits>&
operator>>(std::basic_istream<CharT,Traits>& is, binomial_distribution& bd)
{
bd.read(is);
return is;
}
#endif
/** Returns true if the two distributions will produce the same
sequence of values, given equal generators. */
friend bool operator==(const binomial_distribution& lhs,
const binomial_distribution& rhs)
{
return lhs._t == rhs._t && lhs._p == rhs._p;
}
/** Returns true if the two distributions could produce different
sequences of values, given equal generators. */
friend bool operator!=(const binomial_distribution& lhs,
const binomial_distribution& rhs)
{
return !(lhs == rhs);
}
private:
/// @cond show_private
template<class CharT, class Traits>
void read(std::basic_istream<CharT, Traits>& is) {
param_type parm;
if(is >> parm) {
param(parm);
}
}
bool use_inversion() const
{
// BTRD is safe when np >= 10
return m < 11;
}
// computes the correction factor for the Stirling approximation
// for log(k!)
static RealType fc(IntType k)
{
if(k < 10) return detail::binomial_table<RealType>::table[k];
else {
RealType ikp1 = RealType(1) / (k + 1);
return (RealType(1)/12
- (RealType(1)/360
- (RealType(1)/1260)*(ikp1*ikp1))*(ikp1*ikp1))*ikp1;
}
}
void init()
{
using std::sqrt;
using std::pow;
RealType p = (0.5 < _p)? (1 - _p) : _p;
IntType t = _t;
m = static_cast<IntType>((t+1)*p);
if(use_inversion()) {
q_n = pow((1 - p), static_cast<RealType>(t));
} else {
btrd.r = p/(1-p);
btrd.nr = (t+1)*btrd.r;
btrd.npq = t*p*(1-p);
RealType sqrt_npq = sqrt(btrd.npq);
btrd.b = 1.15 + 2.53 * sqrt_npq;
btrd.a = -0.0873 + 0.0248*btrd.b + 0.01*p;
btrd.c = t*p + 0.5;
btrd.alpha = (2.83 + 5.1/btrd.b) * sqrt_npq;
btrd.v_r = 0.92 - 4.2/btrd.b;
btrd.u_rv_r = 0.86*btrd.v_r;
}
}
template<class URNG>
result_type generate(URNG& urng) const
{
using std::floor;
using std::abs;
using std::log;
while(true) {
RealType u;
RealType v = boost::random::uniform_01<RealType>()(urng);
if(v <= btrd.u_rv_r) {
RealType u = v/btrd.v_r - 0.43;
return static_cast<IntType>(floor(
(2*btrd.a/(0.5 - abs(u)) + btrd.b)*u + btrd.c));
}
if(v >= btrd.v_r) {
u = boost::random::uniform_01<RealType>()(urng) - 0.5;
} else {
u = v/btrd.v_r - 0.93;
u = ((u < 0)? -0.5 : 0.5) - u;
v = boost::random::uniform_01<RealType>()(urng) * btrd.v_r;
}
RealType us = 0.5 - abs(u);
IntType k = static_cast<IntType>(floor((2*btrd.a/us + btrd.b)*u + btrd.c));
if(k < 0 || k > _t) continue;
v = v*btrd.alpha/(btrd.a/(us*us) + btrd.b);
RealType km = abs(k - m);
if(km <= 15) {
RealType f = 1;
if(m < k) {
IntType i = m;
do {
++i;
f = f*(btrd.nr/i - btrd.r);
} while(i != k);
} else if(m > k) {
IntType i = k;
do {
++i;
v = v*(btrd.nr/i - btrd.r);
} while(i != m);
}
if(v <= f) return k;
else continue;
} else {
// final acceptance/rejection
v = log(v);
RealType rho =
(km/btrd.npq)*(((km/3. + 0.625)*km + 1./6)/btrd.npq + 0.5);
RealType t = -km*km/(2*btrd.npq);
if(v < t - rho) return k;
if(v > t + rho) continue;
IntType nm = _t - m + 1;
RealType h = (m + 0.5)*log((m + 1)/(btrd.r*nm))
+ fc(m) + fc(_t - m);
IntType nk = _t - k + 1;
if(v <= h + (_t+1)*log(static_cast<RealType>(nm)/nk)
+ (k + 0.5)*log(nk*btrd.r/(k+1))
- fc(k)
- fc(_t - k))
{
return k;
} else {
continue;
}
}
}
}
template<class URNG>
IntType invert(IntType t, RealType p, URNG& urng) const
{
RealType q = 1 - p;
RealType s = p / q;
RealType a = (t + 1) * s;
RealType r = q_n;
RealType u = boost::random::uniform_01<RealType>()(urng);
IntType x = 0;
while(u > r) {
u = u - r;
++x;
RealType r1 = ((a/x) - s) * r;
// If r gets too small then the round-off error
// becomes a problem. At this point, p(i) is
// decreasing exponentially, so if we just call
// it 0, it's close enough. Note that the
// minimum value of q_n is about 1e-7, so we
// may need to be a little careful to make sure that
// we don't terminate the first time through the loop
// for float. (Hence the test that r is decreasing)
if(r1 < std::numeric_limits<RealType>::epsilon() && r1 < r) {
break;
}
r = r1;
}
return x;
}
// parameters
IntType _t;
RealType _p;
// common data
IntType m;
union {
// for btrd
struct {
RealType r;
RealType nr;
RealType npq;
RealType b;
RealType a;
RealType c;
RealType alpha;
RealType v_r;
RealType u_rv_r;
} btrd;
// for inversion
RealType q_n;
};
/// @endcond
};
} } }
#include <boost/random/detail/enable_warnings.hpp>
#endif
This diff is collapsed.
/* boost random/lognormal_distribution.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Copyright Steven Watanabe 2011
* Distributed under the Boost Software License, Version 1.0. (See
* accompanying file LICENSE_1_0.txt or copy at
* http://www.boost.org/LICENSE_1_0.txt)
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOB_BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP
#define BOB_BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP
#include <boost/config/no_tr1/cmath.hpp> // std::exp, std::sqrt
#include <cassert>
#include <iosfwd>
#include <istream>
#include <boost/limits.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <bob.core/boost/normal_distribution.hpp>
namespace bob {
namespace core {
namespace random {
/**
* Instantiations of class template lognormal_distribution model a
* \random_distribution. Such a distribution produces random numbers
* with \f$\displaystyle p(x) = \frac{1}{x s \sqrt{2\pi}} e^{\frac{-\left(\log(x)-m\right)^2}{2s^2}}\f$
* for x > 0.
*
* @xmlwarning
* This distribution has been updated to match the C++ standard.
* Its behavior has changed from the original
* boost::lognormal_distribution. A backwards compatible
* version is provided in namespace boost.
* @endxmlwarning
*/
template<class RealType = double>
class lognormal_distribution
{
public:
typedef typename bob::core::random::normal_distribution<RealType>::input_type input_type;
typedef RealType result_type;
class param_type
{
public:
typedef lognormal_distribution distribution_type;
/** Constructs the parameters of a lognormal_distribution. */
explicit param_type(RealType m_arg = RealType(0.0),
RealType s_arg = RealType(1.0))
: _m(m_arg), _s(s_arg) {}
/** Returns the "m" parameter of the distribution. */
RealType m() const { return _m; }
/** Returns the "s" parameter of the distribution. */
RealType s() const { return _s; }
/** Writes the parameters to a std::ostream. */
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
{
os << parm._m << " " << parm._s;
return os;
}
/** Reads the parameters from a std::istream. */
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
{
is >> parm._m >> std::ws >> parm._s;
return is;
}
/** Returns true if the two sets of parameters are equal. */
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
{ return lhs._m == rhs._m && lhs._s == rhs._s; }
/** Returns true if the two sets of parameters are different. */
BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
private:
RealType _m;
RealType _s;
};
/**
* Constructs a lognormal_distribution. @c m and @c s are the
* parameters of the distribution.
*/
explicit lognormal_distribution(RealType m_arg = RealType(0.0),
RealType s_arg = RealType(1.0))
: _normal(m_arg, s_arg) {}
/**
* Constructs a lognormal_distribution from its parameters.
*/
explicit lognormal_distribution(const param_type& parm)
: _normal(parm.m(), parm.s()) {}
// compiler-generated copy ctor and assignment operator are fine
/** Returns the m parameter of the distribution. */
RealType m() const { return _normal.mean(); }
/** Returns the s parameter of the distribution. */
RealType s() const { return _normal.sigma(); }
/** Returns the smallest value that the distribution can produce. */
RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
{ return RealType(0); }
/** Returns the largest value that the distribution can produce. */
RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
{ return (std::numeric_limits<RealType>::infinity)(); }
/** Returns the parameters of the distribution. */
param_type param() const { return param_type(m(), s()); }