To check out this repository please hg clone the following URL, or open the URL using EasyMercurial or your preferred Mercurial client.

The primary repository for this project is hosted at https://github.com/sonic-visualiser/sv-dependency-builds .
This repository is a read-only copy which is updated automatically every hour.

Statistics Download as Zip
| Branch: | Tag: | Revision:

root / any / include / boost / math / distributions / chi_squared.hpp @ 160:cff480c41f97

History | View | Annotate | Download (12.3 KB)

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
} // pdf
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