1 /* boost random/non_central_chi_squared_distribution.hpp header file 2 * 3 * Copyright Thijs van den Berg 2014 4 * 5 * Distributed under the Boost Software License, Version 1.0. (See 6 * accompanying file LICENSE_1_0.txt or copy at 7 * http://www.boost.org/LICENSE_1_0.txt) 8 * 9 * See http://www.boost.org for most recent version including documentation. 10 * 11 * $Id$ 12 */ 13 14 #ifndef BOOST_RANDOM_NON_CENTRAL_CHI_SQUARED_DISTRIBUTION_HPP 15 #define BOOST_RANDOM_NON_CENTRAL_CHI_SQUARED_DISTRIBUTION_HPP 16 17 #include <boost/config/no_tr1/cmath.hpp> 18 #include <iosfwd> 19 #include <istream> 20 #include <boost/limits.hpp> 21 #include <boost/random/detail/config.hpp> 22 #include <boost/random/detail/operators.hpp> 23 #include <boost/random/uniform_real_distribution.hpp> 24 #include <boost/random/normal_distribution.hpp> 25 #include <boost/random/chi_squared_distribution.hpp> 26 #include <boost/random/poisson_distribution.hpp> 27 28 namespace boost { 29 namespace random { 30 31 /** 32 * The noncentral chi-squared distribution is a real valued distribution with 33 * two parameter, @c k and @c lambda. The distribution produces values > 0. 34 * 35 * This is the distribution of the sum of squares of k Normal distributed 36 * variates each with variance one and \f$\lambda\f$ the sum of squares of the 37 * normal means. 38 * 39 * The distribution function is 40 * \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$. 41 * where \f$\displaystyle I_\nu(z)\f$ is a modified Bessel function of the 42 * first kind. 43 * 44 * The algorithm is taken from 45 * 46 * @blockquote 47 * "Monte Carlo Methods in Financial Engineering", Paul Glasserman, 48 * 2003, XIII, 596 p, Stochastic Modelling and Applied Probability, Vol. 53, 49 * ISBN 978-0-387-21617-1, p 124, Fig. 3.5. 50 * @endblockquote 51 */ 52 template <typename RealType = double> 53 class non_central_chi_squared_distribution { 54 public: 55 typedef RealType result_type; 56 typedef RealType input_type; 57 58 class param_type { 59 public: 60 typedef non_central_chi_squared_distribution distribution_type; 61 62 /** 63 * Constructs the parameters of a non_central_chi_squared_distribution. 64 * @c k and @c lambda are the parameter of the distribution. 65 * 66 * Requires: k > 0 && lambda > 0 67 */ 68 explicit param_type(RealType k_arg=RealType (1),RealType lambda_arg=RealType (1))69 param_type(RealType k_arg = RealType(1), RealType lambda_arg = RealType(1)) 70 : _k(k_arg), _lambda(lambda_arg) 71 { 72 BOOST_ASSERT(k_arg > RealType(0)); 73 BOOST_ASSERT(lambda_arg > RealType(0)); 74 } 75 76 /** Returns the @c k parameter of the distribution */ k() const77 RealType k() const { return _k; } 78 79 /** Returns the @c lambda parameter of the distribution */ lambda() const80 RealType lambda() const { return _lambda; } 81 82 /** Writes the parameters of the distribution to a @c std::ostream. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,param_type,parm)83 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm) 84 { 85 os << parm._k << ' ' << parm._lambda; 86 return os; 87 } 88 89 /** Reads the parameters of the distribution from a @c std::istream. */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,param_type,parm)90 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm) 91 { 92 is >> parm._k >> std::ws >> parm._lambda; 93 return is; 94 } 95 96 /** Returns true if the parameters have the same values. */ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type,lhs,rhs)97 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs) 98 { return lhs._k == rhs._k && lhs._lambda == rhs._lambda; } 99 100 /** Returns true if the parameters have different values. */ 101 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type) 102 103 private: 104 RealType _k; 105 RealType _lambda; 106 }; 107 108 /** 109 * Construct a @c non_central_chi_squared_distribution object. @c k and 110 * @c lambda are the parameter of the distribution. 111 * 112 * Requires: k > 0 && lambda > 0 113 */ 114 explicit non_central_chi_squared_distribution(RealType k_arg=RealType (1),RealType lambda_arg=RealType (1))115 non_central_chi_squared_distribution(RealType k_arg = RealType(1), RealType lambda_arg = RealType(1)) 116 : _param(k_arg, lambda_arg) 117 { 118 BOOST_ASSERT(k_arg > RealType(0)); 119 BOOST_ASSERT(lambda_arg > RealType(0)); 120 } 121 122 /** 123 * Construct a @c non_central_chi_squared_distribution object from the parameter. 124 */ 125 explicit non_central_chi_squared_distribution(const param_type & parm)126 non_central_chi_squared_distribution(const param_type& parm) 127 : _param( parm ) 128 { } 129 130 /** 131 * Returns a random variate distributed according to the 132 * non central chi squared distribution specified by @c param. 133 */ 134 template<typename URNG> operator ()(URNG & eng,const param_type & parm) const135 RealType operator()(URNG& eng, const param_type& parm) const 136 { return non_central_chi_squared_distribution(parm)(eng); } 137 138 /** 139 * Returns a random variate distributed according to the 140 * non central chi squared distribution. 141 */ 142 template<typename URNG> operator ()(URNG & eng)143 RealType operator()(URNG& eng) 144 { 145 using std::sqrt; 146 if (_param.k() > 1) { 147 boost::random::normal_distribution<RealType> n_dist; 148 boost::random::chi_squared_distribution<RealType> c_dist(_param.k() - RealType(1)); 149 RealType _z = n_dist(eng); 150 RealType _x = c_dist(eng); 151 RealType term1 = _z + sqrt(_param.lambda()); 152 return term1*term1 + _x; 153 } 154 else { 155 boost::random::poisson_distribution<> p_dist(_param.lambda()/RealType(2)); 156 boost::random::poisson_distribution<>::result_type _p = p_dist(eng); 157 boost::random::chi_squared_distribution<RealType> c_dist(_param.k() + RealType(2)*_p); 158 return c_dist(eng); 159 } 160 } 161 162 /** Returns the @c k parameter of the distribution. */ k() const163 RealType k() const { return _param.k(); } 164 165 /** Returns the @c lambda parameter of the distribution. */ lambda() const166 RealType lambda() const { return _param.lambda(); } 167 168 /** Returns the parameters of the distribution. */ param() const169 param_type param() const { return _param; } 170 171 /** Sets parameters of the distribution. */ param(const param_type & parm)172 void param(const param_type& parm) { _param = parm; } 173 174 /** Resets the distribution, so that subsequent uses does not depend on values already produced by it.*/ reset()175 void reset() {} 176 177 /** Returns the smallest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const178 RealType min BOOST_PREVENT_MACRO_SUBSTITUTION() const 179 { return RealType(0); } 180 181 /** Returns the largest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const182 RealType max BOOST_PREVENT_MACRO_SUBSTITUTION() const 183 { return (std::numeric_limits<RealType>::infinity)(); } 184 185 /** Writes the parameters of the distribution to a @c std::ostream. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,non_central_chi_squared_distribution,dist)186 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, non_central_chi_squared_distribution, dist) 187 { 188 os << dist.param(); 189 return os; 190 } 191 192 /** reads the parameters of the distribution from a @c std::istream. */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,non_central_chi_squared_distribution,dist)193 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, non_central_chi_squared_distribution, dist) 194 { 195 param_type parm; 196 if(is >> parm) { 197 dist.param(parm); 198 } 199 return is; 200 } 201 202 /** Returns true if two distributions have the same parameters and produce 203 the same sequence of random numbers given equal generators.*/ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(non_central_chi_squared_distribution,lhs,rhs)204 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(non_central_chi_squared_distribution, lhs, rhs) 205 { return lhs.param() == rhs.param(); } 206 207 /** Returns true if two distributions have different parameters and/or can produce 208 different sequences of random numbers given equal generators.*/ 209 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(non_central_chi_squared_distribution) 210 211 private: 212 213 /// @cond show_private 214 param_type _param; 215 /// @endcond 216 }; 217 218 } // namespace random 219 } // namespace boost 220 221 #endif 222