]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2006, 2007. |
2 | // Copyright Paul A. Bristow 2008, 2010. | |
3 | ||
4 | // Use, modification and distribution are subject to the | |
5 | // Boost Software License, Version 1.0. | |
6 | // (See accompanying file LICENSE_1_0.txt | |
7 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
8 | ||
9 | #ifndef BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP | |
10 | #define BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP | |
11 | ||
12 | #include <boost/math/distributions/fwd.hpp> | |
13 | #include <boost/math/special_functions/gamma.hpp> // for incomplete beta. | |
14 | #include <boost/math/distributions/complement.hpp> // complements | |
15 | #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks | |
16 | #include <boost/math/special_functions/fpclassify.hpp> | |
17 | ||
18 | #include <utility> | |
19 | ||
20 | namespace boost{ namespace math{ | |
21 | ||
22 | template <class RealType = double, class Policy = policies::policy<> > | |
23 | class chi_squared_distribution | |
24 | { | |
25 | public: | |
26 | typedef RealType value_type; | |
27 | typedef Policy policy_type; | |
28 | ||
29 | chi_squared_distribution(RealType i) : m_df(i) | |
30 | { | |
31 | RealType result; | |
32 | detail::check_df( | |
33 | "boost::math::chi_squared_distribution<%1%>::chi_squared_distribution", m_df, &result, Policy()); | |
34 | } // chi_squared_distribution | |
35 | ||
36 | RealType degrees_of_freedom()const | |
37 | { | |
38 | return m_df; | |
39 | } | |
40 | ||
41 | // Parameter estimation: | |
42 | static RealType find_degrees_of_freedom( | |
43 | RealType difference_from_variance, | |
44 | RealType alpha, | |
45 | RealType beta, | |
46 | RealType variance, | |
47 | RealType hint = 100); | |
48 | ||
49 | private: | |
50 | // | |
51 | // Data member: | |
52 | // | |
53 | RealType m_df; // degrees of freedom is a positive real number. | |
54 | }; // class chi_squared_distribution | |
55 | ||
56 | typedef chi_squared_distribution<double> chi_squared; | |
57 | ||
58 | #ifdef BOOST_MSVC | |
59 | #pragma warning(push) | |
60 | #pragma warning(disable:4127) | |
61 | #endif | |
62 | ||
63 | template <class RealType, class Policy> | |
64 | inline const std::pair<RealType, RealType> range(const chi_squared_distribution<RealType, Policy>& /*dist*/) | |
65 | { // Range of permissible values for random variable x. | |
66 | if (std::numeric_limits<RealType>::has_infinity) | |
67 | { | |
68 | return std::pair<RealType, RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()); // 0 to + infinity. | |
69 | } | |
70 | else | |
71 | { | |
72 | using boost::math::tools::max_value; | |
73 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + max. | |
74 | } | |
75 | } | |
76 | ||
77 | #ifdef BOOST_MSVC | |
78 | #pragma warning(pop) | |
79 | #endif | |
80 | ||
81 | template <class RealType, class Policy> | |
82 | inline const std::pair<RealType, RealType> support(const chi_squared_distribution<RealType, Policy>& /*dist*/) | |
83 | { // Range of supported values for random variable x. | |
84 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
85 | return std::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity. | |
86 | } | |
87 | ||
88 | template <class RealType, class Policy> | |
89 | RealType pdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square) | |
90 | { | |
91 | BOOST_MATH_STD_USING // for ADL of std functions | |
92 | RealType degrees_of_freedom = dist.degrees_of_freedom(); | |
93 | // Error check: | |
94 | RealType error_result; | |
95 | ||
96 | static const char* function = "boost::math::pdf(const chi_squared_distribution<%1%>&, %1%)"; | |
97 | ||
98 | if(false == detail::check_df( | |
99 | function, degrees_of_freedom, &error_result, Policy())) | |
100 | return error_result; | |
101 | ||
102 | if((chi_square < 0) || !(boost::math::isfinite)(chi_square)) | |
103 | { | |
104 | return policies::raise_domain_error<RealType>( | |
105 | function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy()); | |
106 | } | |
107 | ||
108 | if(chi_square == 0) | |
109 | { | |
110 | // Handle special cases: | |
111 | if(degrees_of_freedom < 2) | |
112 | { | |
113 | return policies::raise_overflow_error<RealType>( | |
114 | function, 0, Policy()); | |
115 | } | |
116 | else if(degrees_of_freedom == 2) | |
117 | { | |
118 | return 0.5f; | |
119 | } | |
120 | else | |
121 | { | |
122 | return 0; | |
123 | } | |
124 | } | |
125 | ||
126 | return gamma_p_derivative(degrees_of_freedom / 2, chi_square / 2, Policy()) / 2; | |
127 | ||
128 | ||
129 | template <class RealType, class Policy> | |
130 | inline RealType cdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square) | |
131 | { | |
132 | RealType degrees_of_freedom = dist.degrees_of_freedom(); | |
133 | // Error check: | |
134 | RealType error_result; | |
135 | static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)"; | |
136 | ||
137 | if(false == detail::check_df( | |
138 | function, degrees_of_freedom, &error_result, Policy())) | |
139 | return error_result; | |
140 | ||
141 | if((chi_square < 0) || !(boost::math::isfinite)(chi_square)) | |
142 | { | |
143 | return policies::raise_domain_error<RealType>( | |
144 | function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy()); | |
145 | } | |
146 | ||
147 | return boost::math::gamma_p(degrees_of_freedom / 2, chi_square / 2, Policy()); | |
148 | } // cdf | |
149 | ||
150 | template <class RealType, class Policy> | |
151 | inline RealType quantile(const chi_squared_distribution<RealType, Policy>& dist, const RealType& p) | |
152 | { | |
153 | RealType degrees_of_freedom = dist.degrees_of_freedom(); | |
154 | static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)"; | |
155 | // Error check: | |
156 | RealType error_result; | |
157 | if(false == | |
158 | ( | |
159 | detail::check_df(function, degrees_of_freedom, &error_result, Policy()) | |
160 | && detail::check_probability(function, p, &error_result, Policy())) | |
161 | ) | |
162 | return error_result; | |
163 | ||
164 | return 2 * boost::math::gamma_p_inv(degrees_of_freedom / 2, p, Policy()); | |
165 | } // quantile | |
166 | ||
167 | template <class RealType, class Policy> | |
168 | inline RealType cdf(const complemented2_type<chi_squared_distribution<RealType, Policy>, RealType>& c) | |
169 | { | |
170 | RealType const& degrees_of_freedom = c.dist.degrees_of_freedom(); | |
171 | RealType const& chi_square = c.param; | |
172 | static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)"; | |
173 | // Error check: | |
174 | RealType error_result; | |
175 | if(false == detail::check_df( | |
176 | function, degrees_of_freedom, &error_result, Policy())) | |
177 | return error_result; | |
178 | ||
179 | if((chi_square < 0) || !(boost::math::isfinite)(chi_square)) | |
180 | { | |
181 | return policies::raise_domain_error<RealType>( | |
182 | function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy()); | |
183 | } | |
184 | ||
185 | return boost::math::gamma_q(degrees_of_freedom / 2, chi_square / 2, Policy()); | |
186 | } | |
187 | ||
188 | template <class RealType, class Policy> | |
189 | inline RealType quantile(const complemented2_type<chi_squared_distribution<RealType, Policy>, RealType>& c) | |
190 | { | |
191 | RealType const& degrees_of_freedom = c.dist.degrees_of_freedom(); | |
192 | RealType const& q = c.param; | |
193 | static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)"; | |
194 | // Error check: | |
195 | RealType error_result; | |
196 | if(false == ( | |
197 | detail::check_df(function, degrees_of_freedom, &error_result, Policy()) | |
198 | && detail::check_probability(function, q, &error_result, Policy())) | |
199 | ) | |
200 | return error_result; | |
201 | ||
202 | return 2 * boost::math::gamma_q_inv(degrees_of_freedom / 2, q, Policy()); | |
203 | } | |
204 | ||
205 | template <class RealType, class Policy> | |
206 | inline RealType mean(const chi_squared_distribution<RealType, Policy>& dist) | |
207 | { // Mean of Chi-Squared distribution = v. | |
208 | return dist.degrees_of_freedom(); | |
209 | } // mean | |
210 | ||
211 | template <class RealType, class Policy> | |
212 | inline RealType variance(const chi_squared_distribution<RealType, Policy>& dist) | |
213 | { // Variance of Chi-Squared distribution = 2v. | |
214 | return 2 * dist.degrees_of_freedom(); | |
215 | } // variance | |
216 | ||
217 | template <class RealType, class Policy> | |
218 | inline RealType mode(const chi_squared_distribution<RealType, Policy>& dist) | |
219 | { | |
220 | RealType df = dist.degrees_of_freedom(); | |
221 | static const char* function = "boost::math::mode(const chi_squared_distribution<%1%>&)"; | |
222 | // Most sources only define mode for df >= 2, | |
223 | // but for 0 <= df <= 2, the pdf maximum actually occurs at random variate = 0; | |
224 | // So one could extend the definition of mode thus: | |
225 | //if(df < 0) | |
226 | //{ | |
227 | // return policies::raise_domain_error<RealType>( | |
228 | // function, | |
229 | // "Chi-Squared distribution only has a mode for degrees of freedom >= 0, but got degrees of freedom = %1%.", | |
230 | // df, Policy()); | |
231 | //} | |
232 | //return (df <= 2) ? 0 : df - 2; | |
233 | ||
234 | if(df < 2) | |
235 | return policies::raise_domain_error<RealType>( | |
236 | function, | |
237 | "Chi-Squared distribution only has a mode for degrees of freedom >= 2, but got degrees of freedom = %1%.", | |
238 | df, Policy()); | |
239 | return df - 2; | |
240 | } | |
241 | ||
242 | //template <class RealType, class Policy> | |
243 | //inline RealType median(const chi_squared_distribution<RealType, Policy>& dist) | |
244 | //{ // Median is given by Quantile[dist, 1/2] | |
245 | // RealType df = dist.degrees_of_freedom(); | |
246 | // if(df <= 1) | |
247 | // return tools::domain_error<RealType>( | |
248 | // BOOST_CURRENT_FUNCTION, | |
249 | // "The Chi-Squared distribution only has a mode for degrees of freedom >= 2, but got degrees of freedom = %1%.", | |
250 | // df); | |
251 | // return df - RealType(2)/3; | |
252 | //} | |
253 | // Now implemented via quantile(half) in derived accessors. | |
254 | ||
255 | template <class RealType, class Policy> | |
256 | inline RealType skewness(const chi_squared_distribution<RealType, Policy>& dist) | |
257 | { | |
258 | BOOST_MATH_STD_USING // For ADL | |
259 | RealType df = dist.degrees_of_freedom(); | |
260 | return sqrt (8 / df); // == 2 * sqrt(2 / df); | |
261 | } | |
262 | ||
263 | template <class RealType, class Policy> | |
264 | inline RealType kurtosis(const chi_squared_distribution<RealType, Policy>& dist) | |
265 | { | |
266 | RealType df = dist.degrees_of_freedom(); | |
267 | return 3 + 12 / df; | |
268 | } | |
269 | ||
270 | template <class RealType, class Policy> | |
271 | inline RealType kurtosis_excess(const chi_squared_distribution<RealType, Policy>& dist) | |
272 | { | |
273 | RealType df = dist.degrees_of_freedom(); | |
274 | return 12 / df; | |
275 | } | |
276 | ||
277 | // | |
278 | // Parameter estimation comes last: | |
279 | // | |
280 | namespace detail | |
281 | { | |
282 | ||
283 | template <class RealType, class Policy> | |
284 | struct df_estimator | |
285 | { | |
286 | df_estimator(RealType a, RealType b, RealType variance, RealType delta) | |
287 | : alpha(a), beta(b), ratio(delta/variance) | |
288 | { // Constructor | |
289 | } | |
290 | ||
291 | RealType operator()(const RealType& df) | |
292 | { | |
293 | if(df <= tools::min_value<RealType>()) | |
294 | return 1; | |
295 | chi_squared_distribution<RealType, Policy> cs(df); | |
296 | ||
297 | RealType result; | |
298 | if(ratio > 0) | |
299 | { | |
300 | RealType r = 1 + ratio; | |
301 | result = cdf(cs, quantile(complement(cs, alpha)) / r) - beta; | |
302 | } | |
303 | else | |
304 | { // ratio <= 0 | |
305 | RealType r = 1 + ratio; | |
306 | result = cdf(complement(cs, quantile(cs, alpha) / r)) - beta; | |
307 | } | |
308 | return result; | |
309 | } | |
310 | private: | |
311 | RealType alpha; | |
312 | RealType beta; | |
313 | RealType ratio; // Difference from variance / variance, so fractional. | |
314 | }; | |
315 | ||
316 | } // namespace detail | |
317 | ||
318 | template <class RealType, class Policy> | |
319 | RealType chi_squared_distribution<RealType, Policy>::find_degrees_of_freedom( | |
320 | RealType difference_from_variance, | |
321 | RealType alpha, | |
322 | RealType beta, | |
323 | RealType variance, | |
324 | RealType hint) | |
325 | { | |
326 | static const char* function = "boost::math::chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)"; | |
327 | // Check for domain errors: | |
328 | RealType error_result; | |
329 | if(false == | |
330 | detail::check_probability(function, alpha, &error_result, Policy()) | |
331 | && detail::check_probability(function, beta, &error_result, Policy())) | |
332 | { // Either probability is outside 0 to 1. | |
333 | return error_result; | |
334 | } | |
335 | ||
336 | if(hint <= 0) | |
337 | { // No hint given, so guess df = 1. | |
338 | hint = 1; | |
339 | } | |
340 | ||
341 | detail::df_estimator<RealType, Policy> f(alpha, beta, variance, difference_from_variance); | |
342 | tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>()); | |
343 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
344 | std::pair<RealType, RealType> r = | |
345 | tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy()); | |
346 | RealType result = r.first + (r.second - r.first) / 2; | |
347 | if(max_iter >= policies::get_max_root_iterations<Policy>()) | |
348 | { | |
349 | policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" | |
350 | " either there is no answer to how many degrees of freedom are required" | |
351 | " or the answer is infinite. Current best guess is %1%", result, Policy()); | |
352 | } | |
353 | return result; | |
354 | } | |
355 | ||
356 | } // namespace math | |
357 | } // namespace boost | |
358 | ||
359 | // This include must be at the end, *after* the accessors | |
360 | // for this distribution have been defined, in order to | |
361 | // keep compilers that support two-phase lookup happy. | |
362 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
363 | ||
364 | #endif // BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP |