]>
Commit | Line | Data |
---|---|---|
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 | ||
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 |