Boost165GammaDistribution.hpp

Go to the documentation of this file.

/* boost random/gamma_distribution.hpp header file
 *
 * Copyright Jens Maurer 2002
 * 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 BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP
#define BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP

#include <boost/assert.hpp>
#include <boost/config/no_tr1/cmath.hpp>
#include <boost/limits.hpp>
#include <boost/random/detail/config.hpp>
//#include <boost/random/exponential_distribution.hpp> // swapped for chaste copy of exponential v1.65
#include <boost/static_assert.hpp>
#include <iosfwd>
#include <istream>
#include "Boost165ExponentialDistribution.hpp"

namespace boost
{
namespace random
{

    // The algorithm is taken from Knuth

    template <class RealType = double>
    class gamma_distribution_v165
    {
    public:
        typedef RealType input_type;
        typedef RealType result_type;

        class param_type
        {
        public:
            typedef gamma_distribution_v165 distribution_type;

            param_type(const RealType& alpha_arg = RealType(1.0),
                       const RealType& beta_arg = RealType(1.0))
                    : _alpha(alpha_arg), _beta(beta_arg)
            {
            }

            RealType alpha() const { return _alpha; }
            RealType beta() const { return _beta; }

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS

            template <class CharT, class Traits>
            friend std::basic_ostream<CharT, Traits>&
            operator<<(std::basic_ostream<CharT, Traits>& os,
                       const param_type& parm)
            {
                os << parm._alpha << ' ' << parm._beta;
                return os;
            }

            template <class CharT, class Traits>
            friend std::basic_istream<CharT, Traits>&
            operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
            {
                is >> parm._alpha >> std::ws >> parm._beta;
                return is;
            }
#endif

            friend bool operator==(const param_type& lhs, const param_type& rhs)
            {
                return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
            }
            friend bool operator!=(const param_type& lhs, const param_type& rhs)
            {
                return !(lhs == rhs);
            }

        private:
            RealType _alpha;
            RealType _beta;
        };

#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
        BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
#endif

        explicit gamma_distribution_v165(const result_type& alpha_arg = result_type(1.0),
                                         const result_type& beta_arg = result_type(1.0))
                : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
        {
            BOOST_ASSERT(_alpha > result_type(0));
            BOOST_ASSERT(_beta > result_type(0));
            init();
        }

        explicit gamma_distribution_v165(const param_type& parm)
                : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
        {
            init();
        }

        // compiler-generated copy ctor and assignment operator are fine

        RealType alpha() const { return _alpha; }
        RealType beta() const { return _beta; }
        RealType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
        /* Returns the largest value that the distribution can produce. */
        RealType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
        {
            return (std::numeric_limits<RealType>::infinity)();
        }

        param_type param() const { return param_type(_alpha, _beta); }
        void param(const param_type& parm)
        {
            _alpha = parm.alpha();
            _beta = parm.beta();
            init();
        }

        void reset() { _exp.reset(); }

        template <class Engine>
        result_type operator()(Engine& eng)
        {
#ifndef BOOST_NO_STDC_NAMESPACE
            // allow for Koenig lookup
            using std::tan;
            using std::sqrt;
            using std::exp;
            using std::log;
            using std::pow;
#endif
            if (_alpha == result_type(1))
            {
                return _exp(eng) * _beta;
            }
            else if (_alpha > result_type(1))
            {
                // Can we have a boost::mathconst please?
                const result_type pi = result_type(3.14159265358979323846);
                for (;;)
                {
                    result_type y = tan(pi * uniform_01<RealType>()(eng));
                    result_type x = sqrt(result_type(2) * _alpha - result_type(1)) * y
                        + _alpha - result_type(1);
                    if (x <= result_type(0))
                        continue;
                    if (uniform_01<RealType>()(eng) > (result_type(1) + y * y) * exp((_alpha - result_type(1))
                                                                                         * log(x / (_alpha - result_type(1)))
                                                                                     - sqrt(result_type(2) * _alpha
                                                                                            - result_type(1))
                                                                                         * y))
                        continue;
                    return x * _beta;
                }
            }
            else /* alpha < 1.0 */
            {
                for (;;)
                {
                    result_type u = uniform_01<RealType>()(eng);
                    result_type y = _exp(eng);
                    result_type x, q;
                    if (u < _p)
                    {
                        x = exp(-y / _alpha);
                        q = _p * exp(-x);
                    }
                    else
                    {
                        x = result_type(1) + y;
                        q = _p + (result_type(1) - _p) * pow(x, _alpha - result_type(1));
                    }
                    if (u >= q)
                        continue;
                    return x * _beta;
                }
            }
        }

        template <class URNG>
        RealType operator()(URNG& urng, const param_type& parm) const
        {
            return gamma_distribution_v165(parm)(urng);
        }

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS

        template <class CharT, class Traits>
        friend std::basic_ostream<CharT, Traits>&
        operator<<(std::basic_ostream<CharT, Traits>& os,
                   const gamma_distribution_v165& gd)
        {
            os << gd.param();
            return os;
        }

        template <class CharT, class Traits>
        friend std::basic_istream<CharT, Traits>&
        operator>>(std::basic_istream<CharT, Traits>& is, gamma_distribution_v165& gd)
        {
            gd.read(is);
            return is;
        }
#endif

        friend bool operator==(const gamma_distribution_v165& lhs,
                               const gamma_distribution_v165& rhs)
        {
            return lhs._alpha == rhs._alpha
                && lhs._beta == rhs._beta
                && lhs._exp == rhs._exp;
        }

        friend bool operator!=(const gamma_distribution_v165& lhs,
                               const gamma_distribution_v165& rhs)
        {
            return !(lhs == rhs);
        }

    private:

        template <class CharT, class Traits>
        void read(std::basic_istream<CharT, Traits>& is)
        {
            param_type parm;
            if (is >> parm)
            {
                param(parm);
            }
        }

        void init()
        {
#ifndef BOOST_NO_STDC_NAMESPACE
            // allow for Koenig lookup
            using std::exp;
#endif
            _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
        }

        exponential_distribution_v165<RealType> _exp;
        result_type _alpha;
        result_type _beta;
        // some data precomputed from the parameters
        result_type _p;
    };

} // namespace random

} // namespace boost

#endif // BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP