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
20namespace boost{ namespace math{
21
22template <class RealType = double, class Policy = policies::policy<> >
23class chi_squared_distribution
24{
25public:
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
49private:
50 //
51 // Data member:
52 //
53 RealType m_df; // degrees of freedom is a positive real number.
54}; // class chi_squared_distribution
55
56typedef chi_squared_distribution<double> chi_squared;
57
58#ifdef BOOST_MSVC
59#pragma warning(push)
60#pragma warning(disable:4127)
61#endif
62
63template <class RealType, class Policy>
64inline 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
81template <class RealType, class Policy>
82inline 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
88template <class RealType, class Policy>
89RealType 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} // pdf
128
129template <class RealType, class Policy>
130inline 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
150template <class RealType, class Policy>
151inline 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
167template <class RealType, class Policy>
168inline 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
188template <class RealType, class Policy>
189inline 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
205template <class RealType, class Policy>
206inline 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
211template <class RealType, class Policy>
212inline 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
217template <class RealType, class Policy>
218inline 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
255template <class RealType, class Policy>
256inline 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
263template <class RealType, class Policy>
264inline RealType kurtosis(const chi_squared_distribution<RealType, Policy>& dist)
265{
266 RealType df = dist.degrees_of_freedom();
267 return 3 + 12 / df;
268}
269
270template <class RealType, class Policy>
271inline 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//
280namespace detail
281{
282
283template <class RealType, class Policy>
284struct 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 }
310private:
311 RealType alpha;
312 RealType beta;
313 RealType ratio; // Difference from variance / variance, so fractional.
314};
315
316} // namespace detail
317
318template <class RealType, class Policy>
319RealType 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
365

source code of include/boost/math/distributions/chi_squared.hpp