boost/random/non_central_chi_squared_distribution.hpp

0001 /* boost random/non_central_chi_squared_distribution.hpp header file
0002  *
0003  * Copyright Thijs van den Berg 2014
0004  *
0005  * Distributed under the Boost Software License, Version 1.0. (See
0006  * accompanying file LICENSE_1_0.txt or copy at
0007  * http://www.boost.org/LICENSE_1_0.txt)
0008  *
0009  * See http://www.boost.org for most recent version including documentation.
0010  *
0011  * $Id$
0012  */
0013
0014 #ifndef BOOST_RANDOM_NON_CENTRAL_CHI_SQUARED_DISTRIBUTION_HPP
0015 #define BOOST_RANDOM_NON_CENTRAL_CHI_SQUARED_DISTRIBUTION_HPP
0016
0017 #include <boost/config/no_tr1/cmath.hpp>
0018 #include <iosfwd>
0019 #include <istream>
0020 #include <boost/limits.hpp>
0021 #include <boost/random/detail/config.hpp>
0022 #include <boost/random/detail/operators.hpp>
0023 #include <boost/random/uniform_real_distribution.hpp>
0024 #include <boost/random/normal_distribution.hpp>
0025 #include <boost/random/chi_squared_distribution.hpp>
0026 #include <boost/random/poisson_distribution.hpp>
0027
0028 namespace boost {
0029 namespace random {
0030
0031 /**
0032  * The noncentral chi-squared distribution is a real valued distribution with
0033  * two parameter, @c k and @c lambda.  The distribution produces values > 0.
0034  *
0035  * This is the distribution of the sum of squares of k Normal distributed
0036  * variates each with variance one and \f$\lambda\f$ the sum of squares of the
0037  * normal means.
0038  *
0039  * The distribution function is
0040  * \f$\displaystyle P(x) = \frac{1}{2} e^{-(x+\lambda)/2} \left( \frac{x}{\lambda} \right)^{k/4-1/2} I_{k/2-1}( \sqrt{\lambda x} )\f$.
0041  *  where  \f$\displaystyle I_\nu(z)\f$ is a modified Bessel function of the
0042  * first kind.
0043  *
0044  * The algorithm is taken from
0045  *
0046  *  @blockquote
0047  *  "Monte Carlo Methods in Financial Engineering", Paul Glasserman,
0048  *  2003, XIII, 596 p, Stochastic Modelling and Applied Probability, Vol. 53,
0049  *  ISBN 978-0-387-21617-1, p 124, Fig. 3.5.
0050  *  @endblockquote
0051  */
0052 template <typename RealType = double>
0053 class non_central_chi_squared_distribution {
0054 public:
0055     typedef RealType result_type;
0056     typedef RealType input_type;
0057
0058     class param_type {
0059     public:
0060         typedef non_central_chi_squared_distribution distribution_type;
0061
0062         /**
0063          * Constructs the parameters of a non_central_chi_squared_distribution.
0064          * @c k and @c lambda are the parameter of the distribution.
0065          *
0066          * Requires: k > 0 && lambda > 0
0067          */
0068         explicit
0069         param_type(RealType k_arg = RealType(1), RealType lambda_arg = RealType(1))
0070         : _k(k_arg), _lambda(lambda_arg)
0071         {
0072             BOOST_ASSERT(k_arg > RealType(0));
0073             BOOST_ASSERT(lambda_arg > RealType(0));
0074         }
0075
0076         /** Returns the @c k parameter of the distribution */
0077         RealType k() const { return _k; }
0078
0079         /** Returns the @c lambda parameter of the distribution */
0080         RealType lambda() const { return _lambda; }
0081
0082         /** Writes the parameters of the distribution to a @c std::ostream. */
0083         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
0084         {
0085             os << parm._k << ' ' << parm._lambda;
0086             return os;
0087         }
0088
0089         /** Reads the parameters of the distribution from a @c std::istream. */
0090         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
0091         {
0092             is >> parm._k >> std::ws >> parm._lambda;
0093             return is;
0094         }
0095
0096         /** Returns true if the parameters have the same values. */
0097         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
0098         { return lhs._k == rhs._k && lhs._lambda == rhs._lambda; }
0099
0100         /** Returns true if the parameters have different values. */
0101         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
0102
0103     private:
0104         RealType _k;
0105         RealType _lambda;
0106     };
0107
0108     /**
0109      * Construct a @c non_central_chi_squared_distribution object. @c k and
0110      * @c lambda are the parameter of the distribution.
0111      *
0112      * Requires: k > 0 && lambda > 0
0113      */
0114     explicit
0115     non_central_chi_squared_distribution(RealType k_arg = RealType(1), RealType lambda_arg = RealType(1))
0116       : _param(k_arg, lambda_arg)
0117     {
0118         BOOST_ASSERT(k_arg > RealType(0));
0119         BOOST_ASSERT(lambda_arg > RealType(0));
0120     }
0121
0122     /**
0123      * Construct a @c non_central_chi_squared_distribution object from the parameter.
0124      */
0125     explicit
0126     non_central_chi_squared_distribution(const param_type& parm)
0127       : _param( parm )
0128     { }
0129
0130     /**
0131      * Returns a random variate distributed according to the
0132      * non central chi squared distribution specified by @c param.
0133      */
0134     template<typename URNG>
0135     RealType operator()(URNG& eng, const param_type& parm) const
0136     { return non_central_chi_squared_distribution(parm)(eng); }
0137
0138     /**
0139      * Returns a random variate distributed according to the
0140      * non central chi squared distribution.
0141      */
0142     template<typename URNG>
0143     RealType operator()(URNG& eng)
0144     {
0145         using std::sqrt;
0146         if (_param.k() > 1) {
0147             boost::random::normal_distribution<RealType> n_dist;
0148             boost::random::chi_squared_distribution<RealType> c_dist(_param.k() - RealType(1));
0149             RealType _z = n_dist(eng);
0150             RealType _x = c_dist(eng);
0151             RealType term1 = _z + sqrt(_param.lambda());
0152             return term1*term1 + _x;
0153         }
0154         else {
0155             boost::random::poisson_distribution<> p_dist(_param.lambda()/RealType(2));
0156             boost::random::poisson_distribution<>::result_type _p = p_dist(eng);
0157             boost::random::chi_squared_distribution<RealType> c_dist(_param.k() + RealType(2)*_p);
0158             return c_dist(eng);
0159         }
0160     }
0161
0162     /** Returns the @c k parameter of the distribution. */
0163     RealType k() const { return _param.k(); }
0164
0165     /** Returns the @c lambda parameter of the distribution. */
0166     RealType lambda() const { return _param.lambda(); }
0167
0168     /** Returns the parameters of the distribution. */
0169     param_type param() const { return _param; }
0170
0171     /** Sets parameters of the distribution. */
0172     void param(const param_type& parm) { _param = parm; }
0173
0174     /** Resets the distribution, so that subsequent uses does not depend on values already produced by it.*/
0175     void reset() {}
0176
0177     /** Returns the smallest value that the distribution can produce. */
0178     RealType min BOOST_PREVENT_MACRO_SUBSTITUTION() const
0179     { return RealType(0); }
0180
0181     /** Returns the largest value that the distribution can produce. */
0182     RealType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
0183     { return (std::numeric_limits<RealType>::infinity)(); }
0184
0185     /** Writes the parameters of the distribution to a @c std::ostream. */
0186     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, non_central_chi_squared_distribution, dist)
0187     {
0188         os << dist.param();
0189         return os;
0190     }
0191
0192     /** reads the parameters of the distribution from a @c std::istream. */
0193     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, non_central_chi_squared_distribution, dist)
0194     {
0195         param_type parm;
0196         if(is >> parm) {
0197             dist.param(parm);
0198         }
0199         return is;
0200     }
0201
0202     /** Returns true if two distributions have the same parameters and produce
0203         the same sequence of random numbers given equal generators.*/
0204     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(non_central_chi_squared_distribution, lhs, rhs)
0205     { return lhs.param() == rhs.param(); }
0206
0207     /** Returns true if two distributions have different parameters and/or can produce
0208        different sequences of random numbers given equal generators.*/
0209     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(non_central_chi_squared_distribution)
0210
0211 private:
0212
0213     /// @cond show_private
0214     param_type  _param;
0215     /// @endcond
0216 };
0217
0218 } // namespace random
0219 } // namespace boost
0220
0221 #endif