]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2006. |
2 | ||
3 | // Use, modification and distribution are subject to the | |
4 | // Boost Software License, Version 1.0. | |
5 | // (See accompanying file LICENSE_1_0.txt | |
6 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
7 | ||
8 | #ifndef BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP | |
9 | #define BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP | |
10 | ||
11 | #include <boost/math/distributions/fwd.hpp> | |
12 | #include <boost/math/special_functions/beta.hpp> // for incomplete beta. | |
13 | #include <boost/math/distributions/complement.hpp> // complements | |
14 | #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks | |
15 | #include <boost/math/special_functions/fpclassify.hpp> | |
16 | ||
17 | #include <utility> | |
18 | ||
19 | namespace boost{ namespace math{ | |
20 | ||
21 | template <class RealType = double, class Policy = policies::policy<> > | |
22 | class fisher_f_distribution | |
23 | { | |
24 | public: | |
25 | typedef RealType value_type; | |
26 | typedef Policy policy_type; | |
27 | ||
28 | fisher_f_distribution(const RealType& i, const RealType& j) : m_df1(i), m_df2(j) | |
29 | { | |
30 | static const char* function = "fisher_f_distribution<%1%>::fisher_f_distribution"; | |
31 | RealType result; | |
32 | detail::check_df( | |
33 | function, m_df1, &result, Policy()); | |
34 | detail::check_df( | |
35 | function, m_df2, &result, Policy()); | |
36 | } // fisher_f_distribution | |
37 | ||
38 | RealType degrees_of_freedom1()const | |
39 | { | |
40 | return m_df1; | |
41 | } | |
42 | RealType degrees_of_freedom2()const | |
43 | { | |
44 | return m_df2; | |
45 | } | |
46 | ||
47 | private: | |
48 | // | |
49 | // Data members: | |
50 | // | |
51 | RealType m_df1; // degrees of freedom are a real number. | |
52 | RealType m_df2; // degrees of freedom are a real number. | |
53 | }; | |
54 | ||
55 | typedef fisher_f_distribution<double> fisher_f; | |
56 | ||
57 | template <class RealType, class Policy> | |
58 | inline const std::pair<RealType, RealType> range(const fisher_f_distribution<RealType, Policy>& /*dist*/) | |
59 | { // Range of permissible values for random variable x. | |
60 | using boost::math::tools::max_value; | |
61 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); | |
62 | } | |
63 | ||
64 | template <class RealType, class Policy> | |
65 | inline const std::pair<RealType, RealType> support(const fisher_f_distribution<RealType, Policy>& /*dist*/) | |
66 | { // Range of supported values for random variable x. | |
67 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
68 | using boost::math::tools::max_value; | |
69 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); | |
70 | } | |
71 | ||
72 | template <class RealType, class Policy> | |
73 | RealType pdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x) | |
74 | { | |
75 | BOOST_MATH_STD_USING // for ADL of std functions | |
76 | RealType df1 = dist.degrees_of_freedom1(); | |
77 | RealType df2 = dist.degrees_of_freedom2(); | |
78 | // Error check: | |
79 | RealType error_result = 0; | |
80 | static const char* function = "boost::math::pdf(fisher_f_distribution<%1%> const&, %1%)"; | |
81 | if(false == (detail::check_df( | |
82 | function, df1, &error_result, Policy()) | |
83 | && detail::check_df( | |
84 | function, df2, &error_result, Policy()))) | |
85 | return error_result; | |
86 | ||
87 | if((x < 0) || !(boost::math::isfinite)(x)) | |
88 | { | |
89 | return policies::raise_domain_error<RealType>( | |
90 | function, "Random variable parameter was %1%, but must be > 0 !", x, Policy()); | |
91 | } | |
92 | ||
93 | if(x == 0) | |
94 | { | |
95 | // special cases: | |
96 | if(df1 < 2) | |
97 | return policies::raise_overflow_error<RealType>( | |
98 | function, 0, Policy()); | |
99 | else if(df1 == 2) | |
100 | return 1; | |
101 | else | |
102 | return 0; | |
103 | } | |
104 | ||
105 | // | |
106 | // You reach this formula by direct differentiation of the | |
107 | // cdf expressed in terms of the incomplete beta. | |
108 | // | |
109 | // There are two versions so we don't pass a value of z | |
110 | // that is very close to 1 to ibeta_derivative: for some values | |
111 | // of df1 and df2, all the change takes place in this area. | |
112 | // | |
113 | RealType v1x = df1 * x; | |
114 | RealType result; | |
115 | if(v1x > df2) | |
116 | { | |
117 | result = (df2 * df1) / ((df2 + v1x) * (df2 + v1x)); | |
118 | result *= ibeta_derivative(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy()); | |
119 | } | |
120 | else | |
121 | { | |
122 | result = df2 + df1 * x; | |
123 | result = (result * df1 - x * df1 * df1) / (result * result); | |
124 | result *= ibeta_derivative(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy()); | |
125 | } | |
126 | return result; | |
127 | ||
128 | ||
129 | template <class RealType, class Policy> | |
130 | inline RealType cdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x) | |
131 | { | |
132 | static const char* function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)"; | |
133 | RealType df1 = dist.degrees_of_freedom1(); | |
134 | RealType df2 = dist.degrees_of_freedom2(); | |
135 | // Error check: | |
136 | RealType error_result = 0; | |
137 | if(false == detail::check_df( | |
138 | function, df1, &error_result, Policy()) | |
139 | && detail::check_df( | |
140 | function, df2, &error_result, Policy())) | |
141 | return error_result; | |
142 | ||
143 | if((x < 0) || !(boost::math::isfinite)(x)) | |
144 | { | |
145 | return policies::raise_domain_error<RealType>( | |
146 | function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy()); | |
147 | } | |
148 | ||
149 | RealType v1x = df1 * x; | |
150 | // | |
151 | // There are two equivalent formulas used here, the aim is | |
152 | // to prevent the final argument to the incomplete beta | |
153 | // from being too close to 1: for some values of df1 and df2 | |
154 | // the rate of change can be arbitrarily large in this area, | |
155 | // whilst the value we're passing will have lost information | |
156 | // content as a result of being 0.999999something. Better | |
157 | // to switch things around so we're passing 1-z instead. | |
158 | // | |
159 | return v1x > df2 | |
160 | ? boost::math::ibetac(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy()) | |
161 | : boost::math::ibeta(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy()); | |
162 | } // cdf | |
163 | ||
164 | template <class RealType, class Policy> | |
165 | inline RealType quantile(const fisher_f_distribution<RealType, Policy>& dist, const RealType& p) | |
166 | { | |
167 | static const char* function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)"; | |
168 | RealType df1 = dist.degrees_of_freedom1(); | |
169 | RealType df2 = dist.degrees_of_freedom2(); | |
170 | // Error check: | |
171 | RealType error_result = 0; | |
172 | if(false == (detail::check_df( | |
173 | function, df1, &error_result, Policy()) | |
174 | && detail::check_df( | |
175 | function, df2, &error_result, Policy()) | |
176 | && detail::check_probability( | |
177 | function, p, &error_result, Policy()))) | |
178 | return error_result; | |
179 | ||
180 | // With optimizations turned on, gcc wrongly warns about y being used | |
181 | // uninitializated unless we initialize it to something: | |
182 | RealType x, y(0); | |
183 | ||
184 | x = boost::math::ibeta_inv(df1 / 2, df2 / 2, p, &y, Policy()); | |
185 | ||
186 | return df2 * x / (df1 * y); | |
187 | } // quantile | |
188 | ||
189 | template <class RealType, class Policy> | |
190 | inline RealType cdf(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c) | |
191 | { | |
192 | static const char* function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)"; | |
193 | RealType df1 = c.dist.degrees_of_freedom1(); | |
194 | RealType df2 = c.dist.degrees_of_freedom2(); | |
195 | RealType x = c.param; | |
196 | // Error check: | |
197 | RealType error_result = 0; | |
198 | if(false == detail::check_df( | |
199 | function, df1, &error_result, Policy()) | |
200 | && detail::check_df( | |
201 | function, df2, &error_result, Policy())) | |
202 | return error_result; | |
203 | ||
204 | if((x < 0) || !(boost::math::isfinite)(x)) | |
205 | { | |
206 | return policies::raise_domain_error<RealType>( | |
207 | function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy()); | |
208 | } | |
209 | ||
210 | RealType v1x = df1 * x; | |
211 | // | |
212 | // There are two equivalent formulas used here, the aim is | |
213 | // to prevent the final argument to the incomplete beta | |
214 | // from being too close to 1: for some values of df1 and df2 | |
215 | // the rate of change can be arbitrarily large in this area, | |
216 | // whilst the value we're passing will have lost information | |
217 | // content as a result of being 0.999999something. Better | |
218 | // to switch things around so we're passing 1-z instead. | |
219 | // | |
220 | return v1x > df2 | |
221 | ? boost::math::ibeta(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy()) | |
222 | : boost::math::ibetac(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy()); | |
223 | } | |
224 | ||
225 | template <class RealType, class Policy> | |
226 | inline RealType quantile(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c) | |
227 | { | |
228 | static const char* function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)"; | |
229 | RealType df1 = c.dist.degrees_of_freedom1(); | |
230 | RealType df2 = c.dist.degrees_of_freedom2(); | |
231 | RealType p = c.param; | |
232 | // Error check: | |
233 | RealType error_result = 0; | |
234 | if(false == (detail::check_df( | |
235 | function, df1, &error_result, Policy()) | |
236 | && detail::check_df( | |
237 | function, df2, &error_result, Policy()) | |
238 | && detail::check_probability( | |
239 | function, p, &error_result, Policy()))) | |
240 | return error_result; | |
241 | ||
242 | RealType x, y; | |
243 | ||
244 | x = boost::math::ibetac_inv(df1 / 2, df2 / 2, p, &y, Policy()); | |
245 | ||
246 | return df2 * x / (df1 * y); | |
247 | } | |
248 | ||
249 | template <class RealType, class Policy> | |
250 | inline RealType mean(const fisher_f_distribution<RealType, Policy>& dist) | |
251 | { // Mean of F distribution = v. | |
252 | static const char* function = "boost::math::mean(fisher_f_distribution<%1%> const&)"; | |
253 | RealType df1 = dist.degrees_of_freedom1(); | |
254 | RealType df2 = dist.degrees_of_freedom2(); | |
255 | // Error check: | |
256 | RealType error_result = 0; | |
257 | if(false == detail::check_df( | |
258 | function, df1, &error_result, Policy()) | |
259 | && detail::check_df( | |
260 | function, df2, &error_result, Policy())) | |
261 | return error_result; | |
262 | if(df2 <= 2) | |
263 | { | |
264 | return policies::raise_domain_error<RealType>( | |
265 | function, "Second degree of freedom was %1% but must be > 2 in order for the distribution to have a mean.", df2, Policy()); | |
266 | } | |
267 | return df2 / (df2 - 2); | |
268 | } // mean | |
269 | ||
270 | template <class RealType, class Policy> | |
271 | inline RealType variance(const fisher_f_distribution<RealType, Policy>& dist) | |
272 | { // Variance of F distribution. | |
273 | static const char* function = "boost::math::variance(fisher_f_distribution<%1%> const&)"; | |
274 | RealType df1 = dist.degrees_of_freedom1(); | |
275 | RealType df2 = dist.degrees_of_freedom2(); | |
276 | // Error check: | |
277 | RealType error_result = 0; | |
278 | if(false == detail::check_df( | |
279 | function, df1, &error_result, Policy()) | |
280 | && detail::check_df( | |
281 | function, df2, &error_result, Policy())) | |
282 | return error_result; | |
283 | if(df2 <= 4) | |
284 | { | |
285 | return policies::raise_domain_error<RealType>( | |
286 | function, "Second degree of freedom was %1% but must be > 4 in order for the distribution to have a valid variance.", df2, Policy()); | |
287 | } | |
288 | return 2 * df2 * df2 * (df1 + df2 - 2) / (df1 * (df2 - 2) * (df2 - 2) * (df2 - 4)); | |
289 | } // variance | |
290 | ||
291 | template <class RealType, class Policy> | |
292 | inline RealType mode(const fisher_f_distribution<RealType, Policy>& dist) | |
293 | { | |
294 | static const char* function = "boost::math::mode(fisher_f_distribution<%1%> const&)"; | |
295 | RealType df1 = dist.degrees_of_freedom1(); | |
296 | RealType df2 = dist.degrees_of_freedom2(); | |
297 | // Error check: | |
298 | RealType error_result = 0; | |
299 | if(false == detail::check_df( | |
300 | function, df1, &error_result, Policy()) | |
301 | && detail::check_df( | |
302 | function, df2, &error_result, Policy())) | |
303 | return error_result; | |
304 | if(df2 <= 2) | |
305 | { | |
306 | return policies::raise_domain_error<RealType>( | |
307 | function, "Second degree of freedom was %1% but must be > 2 in order for the distribution to have a mode.", df2, Policy()); | |
308 | } | |
309 | return df2 * (df1 - 2) / (df1 * (df2 + 2)); | |
310 | } | |
311 | ||
312 | //template <class RealType, class Policy> | |
313 | //inline RealType median(const fisher_f_distribution<RealType, Policy>& dist) | |
314 | //{ // Median of Fisher F distribution is not defined. | |
315 | // return tools::domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN()); | |
316 | // } // median | |
317 | ||
318 | // Now implemented via quantile(half) in derived accessors. | |
319 | ||
320 | template <class RealType, class Policy> | |
321 | inline RealType skewness(const fisher_f_distribution<RealType, Policy>& dist) | |
322 | { | |
323 | static const char* function = "boost::math::skewness(fisher_f_distribution<%1%> const&)"; | |
324 | BOOST_MATH_STD_USING // ADL of std names | |
325 | // See http://mathworld.wolfram.com/F-Distribution.html | |
326 | RealType df1 = dist.degrees_of_freedom1(); | |
327 | RealType df2 = dist.degrees_of_freedom2(); | |
328 | // Error check: | |
329 | RealType error_result = 0; | |
330 | if(false == detail::check_df( | |
331 | function, df1, &error_result, Policy()) | |
332 | && detail::check_df( | |
333 | function, df2, &error_result, Policy())) | |
334 | return error_result; | |
335 | if(df2 <= 6) | |
336 | { | |
337 | return policies::raise_domain_error<RealType>( | |
338 | function, "Second degree of freedom was %1% but must be > 6 in order for the distribution to have a skewness.", df2, Policy()); | |
339 | } | |
340 | return 2 * (df2 + 2 * df1 - 2) * sqrt((2 * df2 - 8) / (df1 * (df2 + df1 - 2))) / (df2 - 6); | |
341 | } | |
342 | ||
343 | template <class RealType, class Policy> | |
344 | RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist); | |
345 | ||
346 | template <class RealType, class Policy> | |
347 | inline RealType kurtosis(const fisher_f_distribution<RealType, Policy>& dist) | |
348 | { | |
349 | return 3 + kurtosis_excess(dist); | |
350 | } | |
351 | ||
352 | template <class RealType, class Policy> | |
353 | inline RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist) | |
354 | { | |
355 | static const char* function = "boost::math::kurtosis_excess(fisher_f_distribution<%1%> const&)"; | |
356 | // See http://mathworld.wolfram.com/F-Distribution.html | |
357 | RealType df1 = dist.degrees_of_freedom1(); | |
358 | RealType df2 = dist.degrees_of_freedom2(); | |
359 | // Error check: | |
360 | RealType error_result = 0; | |
361 | if(false == detail::check_df( | |
362 | function, df1, &error_result, Policy()) | |
363 | && detail::check_df( | |
364 | function, df2, &error_result, Policy())) | |
365 | return error_result; | |
366 | if(df2 <= 8) | |
367 | { | |
368 | return policies::raise_domain_error<RealType>( | |
369 | function, "Second degree of freedom was %1% but must be > 8 in order for the distribution to have a kutosis.", df2, Policy()); | |
370 | } | |
371 | RealType df2_2 = df2 * df2; | |
372 | RealType df1_2 = df1 * df1; | |
373 | RealType n = -16 + 20 * df2 - 8 * df2_2 + df2_2 * df2 + 44 * df1 - 32 * df2 * df1 + 5 * df2_2 * df1 - 22 * df1_2 + 5 * df2 * df1_2; | |
374 | n *= 12; | |
375 | RealType d = df1 * (df2 - 6) * (df2 - 8) * (df1 + df2 - 2); | |
376 | return n / d; | |
377 | } | |
378 | ||
379 | } // namespace math | |
380 | } // namespace boost | |
381 | ||
382 | // This include must be at the end, *after* the accessors | |
383 | // for this distribution have been defined, in order to | |
384 | // keep compilers that support two-phase lookup happy. | |
385 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
386 | ||
387 | #endif // BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP |