]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |