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