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 / inverse_chi_squared.hpp @ 160:cff480c41f97

History | View | Annotate | Download (15.2 KB)

1
// Copyright John Maddock 2010.
2
// Copyright Paul A. Bristow 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_INVERSE_CHI_SQUARED_HPP
10
#define BOOST_MATH_DISTRIBUTIONS_INVERSE_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> // for complements.
15
#include <boost/math/distributions/detail/common_error_handling.hpp> // for error checks.
16
#include <boost/math/special_functions/fpclassify.hpp> // for isfinite
17

    
18
// See http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution
19
// for definitions of this scaled version.
20
// See http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
21
// for unscaled version.
22

    
23
// http://reference.wolfram.com/mathematica/ref/InverseChiSquareDistribution.html
24
// Weisstein, Eric W. "Inverse Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource.
25
// http://mathworld.wolfram.com/InverseChi-SquaredDistribution.html
26

    
27
#include <utility>
28

    
29
namespace boost{ namespace math{
30

    
31
namespace detail
32
{
33
  template <class RealType, class Policy>
34
  inline bool check_inverse_chi_squared( // Check both distribution parameters.
35
        const char* function,
36
        RealType degrees_of_freedom, // degrees_of_freedom (aka nu).
37
        RealType scale,  // scale (aka sigma^2)
38
        RealType* result,
39
        const Policy& pol)
40
  {
41
     return check_scale(function, scale, result, pol)
42
       && check_df(function, degrees_of_freedom,
43
       result, pol);
44
  } // bool check_inverse_chi_squared
45
} // namespace detail
46

    
47
template <class RealType = double, class Policy = policies::policy<> >
48
class inverse_chi_squared_distribution
49
{
50
public:
51
   typedef RealType value_type;
52
   typedef Policy policy_type;
53

    
54
   inverse_chi_squared_distribution(RealType df, RealType l_scale) : m_df(df), m_scale (l_scale)
55
   {
56
      RealType result;
57
      detail::check_df(
58
         "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
59
         m_df, &result, Policy())
60
         && detail::check_scale(
61
"boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
62
         m_scale, &result,  Policy());
63
   } // inverse_chi_squared_distribution constructor 
64

    
65
   inverse_chi_squared_distribution(RealType df = 1) : m_df(df)
66
   {
67
      RealType result;
68
      m_scale = 1 / m_df ; // Default scale = 1 / degrees of freedom (Wikipedia definition 1).
69
      detail::check_df(
70
         "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
71
         m_df, &result, Policy());
72
   } // inverse_chi_squared_distribution
73

    
74
   RealType degrees_of_freedom()const
75
   {
76
      return m_df; // aka nu
77
   }
78
   RealType scale()const
79
   {
80
      return m_scale;  // aka xi
81
   }
82

    
83
   // Parameter estimation:  NOT implemented yet.
84
   //static RealType find_degrees_of_freedom(
85
   //   RealType difference_from_variance,
86
   //   RealType alpha,
87
   //   RealType beta,
88
   //   RealType variance,
89
   //   RealType hint = 100);
90

    
91
private:
92
   // Data members:
93
   RealType m_df;  // degrees of freedom are treated as a real number.
94
   RealType m_scale;  // distribution scale.
95

    
96
}; // class chi_squared_distribution
97

    
98
typedef inverse_chi_squared_distribution<double> inverse_chi_squared;
99

    
100
template <class RealType, class Policy>
101
inline const std::pair<RealType, RealType> range(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
102
{  // Range of permissible values for random variable x.
103
   using boost::math::tools::max_value;
104
   return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + infinity.
105
}
106

    
107
template <class RealType, class Policy>
108
inline const std::pair<RealType, RealType> support(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
109
{  // Range of supported values for random variable x.
110
   // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
111
   return std::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity.
112
}
113

    
114
template <class RealType, class Policy>
115
RealType pdf(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
116
{
117
   BOOST_MATH_STD_USING  // for ADL of std functions.
118
   RealType df = dist.degrees_of_freedom();
119
   RealType scale = dist.scale();
120
   RealType error_result;
121

    
122
   static const char* function = "boost::math::pdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
123

    
124
   if(false == detail::check_inverse_chi_squared
125
     (function, df, scale, &error_result, Policy())
126
     )
127
   { // Bad distribution.
128
      return error_result;
129
   }
130
   if((x < 0) || !(boost::math::isfinite)(x))
131
   { // Bad x.
132
      return policies::raise_domain_error<RealType>(
133
         function, "inverse Chi Square parameter was %1%, but must be >= 0 !", x, Policy());
134
   }
135

    
136
   if(x == 0)
137
   { // Treat as special case.
138
     return 0;
139
   }
140
   // Wikipedia scaled inverse chi sq (df, scale) related to inv gamma (df/2, df * scale /2) 
141
   // so use inverse gamma pdf with shape = df/2, scale df * scale /2 
142
   // RealType shape = df /2; // inv_gamma shape
143
   // RealType scale = df * scale/2; // inv_gamma scale
144
   // RealType result = gamma_p_derivative(shape, scale / x, Policy()) * scale / (x * x);
145
   RealType result = df * scale/2 / x;
146
   if(result < tools::min_value<RealType>())
147
      return 0; // Random variable is near enough infinite.
148
   result = gamma_p_derivative(df/2, result, Policy()) * df * scale/2;
149
   if(result != 0) // prevent 0 / 0,  gamma_p_derivative -> 0 faster than x^2
150
      result /= (x * x);
151
   return result;
152
} // pdf
153

    
154
template <class RealType, class Policy>
155
inline RealType cdf(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
156
{
157
   static const char* function = "boost::math::cdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
158
   RealType df = dist.degrees_of_freedom();
159
   RealType scale = dist.scale();
160
   RealType error_result;
161

    
162
   if(false ==
163
       detail::check_inverse_chi_squared(function, df, scale, &error_result, Policy())
164
     )
165
   { // Bad distribution.
166
      return error_result;
167
   }
168
   if((x < 0) || !(boost::math::isfinite)(x))
169
   { // Bad x.
170
      return policies::raise_domain_error<RealType>(
171
         function, "inverse Chi Square parameter was %1%, but must be >= 0 !", x, Policy());
172
   }
173
   if (x == 0)
174
   { // Treat zero as a special case.
175
     return 0;
176
   }
177
   // RealType shape = df /2; // inv_gamma shape,
178
   // RealType scale = df * scale/2; // inv_gamma scale,
179
   // result = boost::math::gamma_q(shape, scale / x, Policy()); // inverse_gamma code.
180
   return boost::math::gamma_q(df / 2, (df * (scale / 2)) / x, Policy());
181
} // cdf
182

    
183
template <class RealType, class Policy>
184
inline RealType quantile(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
185
{
186
   using boost::math::gamma_q_inv;
187
   RealType df = dist.degrees_of_freedom();
188
   RealType scale = dist.scale();
189

    
190
   static const char* function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
191
   // Error check:
192
   RealType error_result;
193
   if(false == detail::check_df(
194
         function, df, &error_result, Policy())
195
         && detail::check_probability(
196
            function, p, &error_result, Policy()))
197
   {
198
      return error_result;
199
   }
200
   if(false == detail::check_probability(
201
            function, p, &error_result, Policy()))
202
   {
203
      return error_result;
204
   }
205
   // RealType shape = df /2; // inv_gamma shape,
206
   // RealType scale = df * scale/2; // inv_gamma scale,
207
   // result = scale / gamma_q_inv(shape, p, Policy());
208
      RealType result = gamma_q_inv(df /2, p, Policy());
209
      if(result == 0)
210
         return policies::raise_overflow_error<RealType, Policy>(function, "Random variable is infinite.", Policy());
211
      result = df * (scale / 2) / result;
212
      return result;
213
} // quantile
214

    
215
template <class RealType, class Policy>
216
inline RealType cdf(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
217
{
218
   using boost::math::gamma_q_inv;
219
   RealType const& df = c.dist.degrees_of_freedom();
220
   RealType const& scale = c.dist.scale();
221
   RealType const& x = c.param;
222
   static const char* function = "boost::math::cdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
223
   // Error check:
224
   RealType error_result;
225
   if(false == detail::check_df(
226
         function, df, &error_result, Policy()))
227
   {
228
      return error_result;
229
   }
230
   if (x == 0)
231
   { // Treat zero as a special case.
232
     return 1;
233
   }
234
   if((x < 0) || !(boost::math::isfinite)(x))
235
   {
236
      return policies::raise_domain_error<RealType>(
237
         function, "inverse Chi Square parameter was %1%, but must be > 0 !", x, Policy());
238
   }
239
   // RealType shape = df /2; // inv_gamma shape,
240
   // RealType scale = df * scale/2; // inv_gamma scale,
241
   // result = gamma_p(shape, scale/c.param, Policy()); use inv_gamma.
242

    
243
   return gamma_p(df / 2, (df * scale/2) / x, Policy()); // OK
244
} // cdf(complemented
245

    
246
template <class RealType, class Policy>
247
inline RealType quantile(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
248
{
249
   using boost::math::gamma_q_inv;
250

    
251
   RealType const& df = c.dist.degrees_of_freedom();
252
   RealType const& scale = c.dist.scale();
253
   RealType const& q = c.param;
254
   static const char* function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
255
   // Error check:
256
   RealType error_result;
257
   if(false == detail::check_df(function, df, &error_result, Policy()))
258
   {
259
      return error_result;
260
   }
261
   if(false == detail::check_probability(function, q, &error_result, Policy()))
262
   {
263
      return error_result;
264
   }
265
   // RealType shape = df /2; // inv_gamma shape,
266
   // RealType scale = df * scale/2; // inv_gamma scale,
267
   // result = scale / gamma_p_inv(shape, q, Policy());  // using inv_gamma.
268
   RealType result = gamma_p_inv(df/2, q, Policy());
269
   if(result == 0)
270
      return policies::raise_overflow_error<RealType, Policy>(function, "Random variable is infinite.", Policy());
271
   result = (df * scale / 2) / result;
272
   return result;
273
} // quantile(const complement
274

    
275
template <class RealType, class Policy>
276
inline RealType mean(const inverse_chi_squared_distribution<RealType, Policy>& dist)
277
{ // Mean of inverse Chi-Squared distribution.
278
   RealType df = dist.degrees_of_freedom();
279
   RealType scale = dist.scale();
280

    
281
   static const char* function = "boost::math::mean(const inverse_chi_squared_distribution<%1%>&)";
282
   if(df <= 2)
283
      return policies::raise_domain_error<RealType>(
284
         function,
285
         "inverse Chi-Squared distribution only has a mode for degrees of freedom > 2, but got degrees of freedom = %1%.",
286
         df, Policy());
287
  return (df * scale) / (df - 2);
288
} // mean
289

    
290
template <class RealType, class Policy>
291
inline RealType variance(const inverse_chi_squared_distribution<RealType, Policy>& dist)
292
{ // Variance of inverse Chi-Squared distribution.
293
   RealType df = dist.degrees_of_freedom();
294
   RealType scale = dist.scale();
295
   static const char* function = "boost::math::variance(const inverse_chi_squared_distribution<%1%>&)";
296
   if(df <= 4)
297
   {
298
      return policies::raise_domain_error<RealType>(
299
         function,
300
         "inverse Chi-Squared distribution only has a variance for degrees of freedom > 4, but got degrees of freedom = %1%.",
301
         df, Policy());
302
   }
303
   return 2 * df * df * scale * scale / ((df - 2)*(df - 2) * (df - 4));
304
} // variance
305

    
306
template <class RealType, class Policy>
307
inline RealType mode(const inverse_chi_squared_distribution<RealType, Policy>& dist)
308
{ // mode is not defined in Mathematica.
309
  // See Discussion section http://en.wikipedia.org/wiki/Talk:Scaled-inverse-chi-square_distribution
310
  // for origin of the formula used below.
311

    
312
   RealType df = dist.degrees_of_freedom();
313
   RealType scale = dist.scale();
314
   static const char* function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
315
   if(df < 0)
316
      return policies::raise_domain_error<RealType>(
317
         function,
318
         "inverse Chi-Squared distribution only has a mode for degrees of freedom >= 0, but got degrees of freedom = %1%.",
319
         df, Policy());
320
   return (df * scale) / (df + 2);
321
}
322

    
323
//template <class RealType, class Policy>
324
//inline RealType median(const inverse_chi_squared_distribution<RealType, Policy>& dist)
325
//{ // Median is given by Quantile[dist, 1/2]
326
//   RealType df = dist.degrees_of_freedom();
327
//   if(df <= 1)
328
//      return tools::domain_error<RealType>(
329
//         BOOST_CURRENT_FUNCTION,
330
//         "The inverse_Chi-Squared distribution only has a median for degrees of freedom >= 0, but got degrees of freedom = %1%.",
331
//         df);
332
//   return df;
333
//}
334
// Now implemented via quantile(half) in derived accessors.
335

    
336
template <class RealType, class Policy>
337
inline RealType skewness(const inverse_chi_squared_distribution<RealType, Policy>& dist)
338
{
339
   BOOST_MATH_STD_USING // For ADL
340
   RealType df = dist.degrees_of_freedom();
341
   static const char* function = "boost::math::skewness(const inverse_chi_squared_distribution<%1%>&)";
342
   if(df <= 6)
343
      return policies::raise_domain_error<RealType>(
344
         function,
345
         "inverse Chi-Squared distribution only has a skewness for degrees of freedom > 6, but got degrees of freedom = %1%.",
346
         df, Policy());
347

    
348
   return 4 * sqrt (2 * (df - 4)) / (df - 6);  // Not a function of scale.
349
}
350

    
351
template <class RealType, class Policy>
352
inline RealType kurtosis(const inverse_chi_squared_distribution<RealType, Policy>& dist)
353
{
354
   RealType df = dist.degrees_of_freedom();
355
   static const char* function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
356
   if(df <= 8)
357
      return policies::raise_domain_error<RealType>(
358
         function,
359
         "inverse Chi-Squared distribution only has a kurtosis for degrees of freedom > 8, but got degrees of freedom = %1%.",
360
         df, Policy());
361

    
362
   return kurtosis_excess(dist) + 3;
363
}
364

    
365
template <class RealType, class Policy>
366
inline RealType kurtosis_excess(const inverse_chi_squared_distribution<RealType, Policy>& dist)
367
{
368
   RealType df = dist.degrees_of_freedom();
369
   static const char* function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
370
   if(df <= 8)
371
      return policies::raise_domain_error<RealType>(
372
         function,
373
         "inverse Chi-Squared distribution only has a kurtosis excess for degrees of freedom > 8, but got degrees of freedom = %1%.",
374
         df, Policy());
375

    
376
   return 12 * (5 * df - 22) / ((df - 6 )*(df - 8));  // Not a function of scale.
377
}
378

    
379
//
380
// Parameter estimation comes last:
381
//
382

    
383
} // namespace math
384
} // namespace boost
385

    
386
// This include must be at the end, *after* the accessors
387
// for this distribution have been defined, in order to
388
// keep compilers that support two-phase lookup happy.
389
#include <boost/math/distributions/detail/derived_accessors.hpp>
390

    
391
#endif // BOOST_MATH_DISTRIBUTIONS_INVERSE_CHI_SQUARED_HPP