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 using value_type = RealType;
27 using policy_type = Policy;
28
29 explicit 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
56using chi_squared = chi_squared_distribution<double>;
57
58#ifdef __cpp_deduction_guides
59template <class RealType>
60chi_squared_distribution(RealType)->chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>;
61#endif
62
63#ifdef _MSC_VER
64#pragma warning(push)
65#pragma warning(disable:4127)
66#endif
67
68template <class RealType, class Policy>
69inline std::pair<RealType, RealType> range(const chi_squared_distribution<RealType, Policy>& /*dist*/)
70{ // Range of permissible values for random variable x.
71 if (std::numeric_limits<RealType>::has_infinity)
72 {
73 return std::pair<RealType, RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()); // 0 to + infinity.
74 }
75 else
76 {
77 using boost::math::tools::max_value;
78 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + max.
79 }
80}
81
82#ifdef _MSC_VER
83#pragma warning(pop)
84#endif
85
86template <class RealType, class Policy>
87inline std::pair<RealType, RealType> support(const chi_squared_distribution<RealType, Policy>& /*dist*/)
88{ // Range of supported values for random variable x.
89 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
90 return std::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity.
91}
92
93template <class RealType, class Policy>
94RealType pdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square)
95{
96 BOOST_MATH_STD_USING // for ADL of std functions
97 RealType degrees_of_freedom = dist.degrees_of_freedom();
98 // Error check:
99 RealType error_result;
100
101 static const char* function = "boost::math::pdf(const chi_squared_distribution<%1%>&, %1%)";
102
103 if(false == detail::check_df(
104 function, degrees_of_freedom, &error_result, Policy()))
105 return error_result;
106
107 if((chi_square < 0) || !(boost::math::isfinite)(chi_square))
108 {
109 return policies::raise_domain_error<RealType>(
110 function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy());
111 }
112
113 if(chi_square == 0)
114 {
115 // Handle special cases:
116 if(degrees_of_freedom < 2)
117 {
118 return policies::raise_overflow_error<RealType>(
119 function, 0, Policy());
120 }
121 else if(degrees_of_freedom == 2)
122 {
123 return 0.5f;
124 }
125 else
126 {
127 return 0;
128 }
129 }
130
131 return gamma_p_derivative(degrees_of_freedom / 2, chi_square / 2, Policy()) / 2;
132} // pdf
133
134template <class RealType, class Policy>
135inline RealType cdf(const chi_squared_distribution<RealType, Policy>& dist, const RealType& chi_square)
136{
137 RealType degrees_of_freedom = dist.degrees_of_freedom();
138 // Error check:
139 RealType error_result;
140 static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)";
141
142 if(false == detail::check_df(
143 function, degrees_of_freedom, &error_result, Policy()))
144 return error_result;
145
146 if((chi_square < 0) || !(boost::math::isfinite)(chi_square))
147 {
148 return policies::raise_domain_error<RealType>(
149 function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy());
150 }
151
152 return boost::math::gamma_p(degrees_of_freedom / 2, chi_square / 2, Policy());
153} // cdf
154
155template <class RealType, class Policy>
156inline RealType quantile(const chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
157{
158 RealType degrees_of_freedom = dist.degrees_of_freedom();
159 static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)";
160 // Error check:
161 RealType error_result;
162 if(false ==
163 (
164 detail::check_df(function, degrees_of_freedom, &error_result, Policy())
165 && detail::check_probability(function, p, &error_result, Policy()))
166 )
167 return error_result;
168
169 return 2 * boost::math::gamma_p_inv(degrees_of_freedom / 2, p, Policy());
170} // quantile
171
172template <class RealType, class Policy>
173inline RealType cdf(const complemented2_type<chi_squared_distribution<RealType, Policy>, RealType>& c)
174{
175 RealType const& degrees_of_freedom = c.dist.degrees_of_freedom();
176 RealType const& chi_square = c.param;
177 static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)";
178 // Error check:
179 RealType error_result;
180 if(false == detail::check_df(
181 function, degrees_of_freedom, &error_result, Policy()))
182 return error_result;
183
184 if((chi_square < 0) || !(boost::math::isfinite)(chi_square))
185 {
186 return policies::raise_domain_error<RealType>(
187 function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy());
188 }
189
190 return boost::math::gamma_q(degrees_of_freedom / 2, chi_square / 2, Policy());
191}
192
193template <class RealType, class Policy>
194inline RealType quantile(const complemented2_type<chi_squared_distribution<RealType, Policy>, RealType>& c)
195{
196 RealType const& degrees_of_freedom = c.dist.degrees_of_freedom();
197 RealType const& q = c.param;
198 static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)";
199 // Error check:
200 RealType error_result;
201 if(false == (
202 detail::check_df(function, degrees_of_freedom, &error_result, Policy())
203 && detail::check_probability(function, q, &error_result, Policy()))
204 )
205 return error_result;
206
207 return 2 * boost::math::gamma_q_inv(degrees_of_freedom / 2, q, Policy());
208}
209
210template <class RealType, class Policy>
211inline RealType mean(const chi_squared_distribution<RealType, Policy>& dist)
212{ // Mean of Chi-Squared distribution = v.
213 return dist.degrees_of_freedom();
214} // mean
215
216template <class RealType, class Policy>
217inline RealType variance(const chi_squared_distribution<RealType, Policy>& dist)
218{ // Variance of Chi-Squared distribution = 2v.
219 return 2 * dist.degrees_of_freedom();
220} // variance
221
222template <class RealType, class Policy>
223inline RealType mode(const chi_squared_distribution<RealType, Policy>& dist)
224{
225 RealType df = dist.degrees_of_freedom();
226 static const char* function = "boost::math::mode(const chi_squared_distribution<%1%>&)";
227
228 if(df < 2)
229 return policies::raise_domain_error<RealType>(
230 function,
231 "Chi-Squared distribution only has a mode for degrees of freedom >= 2, but got degrees of freedom = %1%.",
232 df, Policy());
233 return df - 2;
234}
235
236template <class RealType, class Policy>
237inline RealType skewness(const chi_squared_distribution<RealType, Policy>& dist)
238{
239 BOOST_MATH_STD_USING // For ADL
240 RealType df = dist.degrees_of_freedom();
241 return sqrt (8 / df);
242}
243
244template <class RealType, class Policy>
245inline RealType kurtosis(const chi_squared_distribution<RealType, Policy>& dist)
246{
247 RealType df = dist.degrees_of_freedom();
248 return 3 + 12 / df;
249}
250
251template <class RealType, class Policy>
252inline RealType kurtosis_excess(const chi_squared_distribution<RealType, Policy>& dist)
253{
254 RealType df = dist.degrees_of_freedom();
255 return 12 / df;
256}
257
258//
259// Parameter estimation comes last:
260//
261namespace detail
262{
263
264template <class RealType, class Policy>
265struct df_estimator
266{
267 df_estimator(RealType a, RealType b, RealType variance, RealType delta)
268 : alpha(a), beta(b), ratio(delta/variance)
269 { // Constructor
270 }
271
272 RealType operator()(const RealType& df)
273 {
274 if(df <= tools::min_value<RealType>())
275 return 1;
276 chi_squared_distribution<RealType, Policy> cs(df);
277
278 RealType result;
279 if(ratio > 0)
280 {
281 RealType r = 1 + ratio;
282 result = cdf(cs, quantile(complement(cs, alpha)) / r) - beta;
283 }
284 else
285 { // ratio <= 0
286 RealType r = 1 + ratio;
287 result = cdf(complement(cs, quantile(cs, alpha) / r)) - beta;
288 }
289 return result;
290 }
291private:
292 RealType alpha;
293 RealType beta;
294 RealType ratio; // Difference from variance / variance, so fractional.
295};
296
297} // namespace detail
298
299template <class RealType, class Policy>
300RealType chi_squared_distribution<RealType, Policy>::find_degrees_of_freedom(
301 RealType difference_from_variance,
302 RealType alpha,
303 RealType beta,
304 RealType variance,
305 RealType hint)
306{
307 static const char* function = "boost::math::chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)";
308 // Check for domain errors:
309 RealType error_result;
310 if(false ==
311 detail::check_probability(function, alpha, &error_result, Policy())
312 && detail::check_probability(function, beta, &error_result, Policy()))
313 { // Either probability is outside 0 to 1.
314 return error_result;
315 }
316
317 if(hint <= 0)
318 { // No hint given, so guess df = 1.
319 hint = 1;
320 }
321
322 detail::df_estimator<RealType, Policy> f(alpha, beta, variance, difference_from_variance);
323 tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
324 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
325 std::pair<RealType, RealType> r =
326 tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
327 RealType result = r.first + (r.second - r.first) / 2;
328 if(max_iter >= policies::get_max_root_iterations<Policy>())
329 {
330 policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
331 " either there is no answer to how many degrees of freedom are required"
332 " or the answer is infinite. Current best guess is %1%", result, Policy());
333 }
334 return result;
335}
336
337} // namespace math
338} // namespace boost
339
340// This include must be at the end, *after* the accessors
341// for this distribution have been defined, in order to
342// keep compilers that support two-phase lookup happy.
343#include <boost/math/distributions/detail/derived_accessors.hpp>
344
345#endif // BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP
346

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