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 |
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 */ |
77 | RealType k() const { return _k; } |
78 | |
79 | /** Returns the @c lambda parameter of the distribution */ |
80 | RealType lambda() const { return _lambda; } |
81 | |
82 | /** Writes the parameters of the distribution to a @c std::ostream. */ |
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. */ |
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. */ |
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 |
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 |
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> |
135 | 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> |
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. */ |
163 | RealType k() const { return _param.k(); } |
164 | |
165 | /** Returns the @c lambda parameter of the distribution. */ |
166 | RealType lambda() const { return _param.lambda(); } |
167 | |
168 | /** Returns the parameters of the distribution. */ |
169 | param_type param() const { return _param; } |
170 | |
171 | /** Sets parameters of the distribution. */ |
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.*/ |
175 | void reset() {} |
176 | |
177 | /** Returns the smallest value that the distribution can produce. */ |
178 | RealType min BOOST_PREVENT_MACRO_SUBSTITUTION() const |
179 | { return RealType(0); } |
180 | |
181 | /** Returns the largest value that the distribution can produce. */ |
182 | 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. */ |
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. */ |
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.*/ |
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 | |