]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // boost\math\distributions\non_central_f.hpp |
2 | ||
3 | // Copyright John Maddock 2008. | |
4 | ||
5 | // Use, modification and distribution are subject to the | |
6 | // Boost Software License, Version 1.0. | |
7 | // (See accompanying file LICENSE_1_0.txt | |
8 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
9 | ||
10 | #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP | |
11 | #define BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP | |
12 | ||
13 | #include <boost/math/distributions/non_central_beta.hpp> | |
14 | #include <boost/math/distributions/detail/generic_mode.hpp> | |
15 | #include <boost/math/special_functions/pow.hpp> | |
16 | ||
17 | namespace boost | |
18 | { | |
19 | namespace math | |
20 | { | |
21 | template <class RealType = double, class Policy = policies::policy<> > | |
22 | class non_central_f_distribution | |
23 | { | |
24 | public: | |
25 | typedef RealType value_type; | |
26 | typedef Policy policy_type; | |
27 | ||
28 | non_central_f_distribution(RealType v1_, RealType v2_, RealType lambda) : v1(v1_), v2(v2_), ncp(lambda) | |
29 | { | |
30 | const char* function = "boost::math::non_central_f_distribution<%1%>::non_central_f_distribution(%1%,%1%)"; | |
31 | RealType r; | |
32 | detail::check_df( | |
33 | function, | |
34 | v1, &r, Policy()); | |
35 | detail::check_df( | |
36 | function, | |
37 | v2, &r, Policy()); | |
38 | detail::check_non_centrality( | |
39 | function, | |
40 | lambda, | |
41 | &r, | |
42 | Policy()); | |
43 | } // non_central_f_distribution constructor. | |
44 | ||
45 | RealType degrees_of_freedom1()const | |
46 | { | |
47 | return v1; | |
48 | } | |
49 | RealType degrees_of_freedom2()const | |
50 | { | |
51 | return v2; | |
52 | } | |
53 | RealType non_centrality() const | |
54 | { // Private data getter function. | |
55 | return ncp; | |
56 | } | |
57 | private: | |
58 | // Data member, initialized by constructor. | |
59 | RealType v1; // alpha. | |
60 | RealType v2; // beta. | |
61 | RealType ncp; // non-centrality parameter | |
62 | }; // template <class RealType, class Policy> class non_central_f_distribution | |
63 | ||
64 | typedef non_central_f_distribution<double> non_central_f; // Reserved name of type double. | |
65 | ||
66 | // Non-member functions to give properties of the distribution. | |
67 | ||
68 | template <class RealType, class Policy> | |
69 | inline const std::pair<RealType, RealType> range(const non_central_f_distribution<RealType, Policy>& /* dist */) | |
70 | { // Range of permissible values for random variable k. | |
71 | using boost::math::tools::max_value; | |
72 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); | |
73 | } | |
74 | ||
75 | template <class RealType, class Policy> | |
76 | inline const std::pair<RealType, RealType> support(const non_central_f_distribution<RealType, Policy>& /* dist */) | |
77 | { // Range of supported values for random variable k. | |
78 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
79 | using boost::math::tools::max_value; | |
80 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); | |
81 | } | |
82 | ||
83 | template <class RealType, class Policy> | |
84 | inline RealType mean(const non_central_f_distribution<RealType, Policy>& dist) | |
85 | { | |
86 | const char* function = "mean(non_central_f_distribution<%1%> const&)"; | |
87 | RealType v1 = dist.degrees_of_freedom1(); | |
88 | RealType v2 = dist.degrees_of_freedom2(); | |
89 | RealType l = dist.non_centrality(); | |
90 | RealType r; | |
91 | if(!detail::check_df( | |
92 | function, | |
93 | v1, &r, Policy()) | |
94 | || | |
95 | !detail::check_df( | |
96 | function, | |
97 | v2, &r, Policy()) | |
98 | || | |
99 | !detail::check_non_centrality( | |
100 | function, | |
101 | l, | |
102 | &r, | |
103 | Policy())) | |
104 | return r; | |
105 | if(v2 <= 2) | |
106 | return policies::raise_domain_error( | |
107 | function, | |
108 | "Second degrees of freedom parameter was %1%, but must be > 2 !", | |
109 | v2, Policy()); | |
110 | return v2 * (v1 + l) / (v1 * (v2 - 2)); | |
111 | } // mean | |
112 | ||
113 | template <class RealType, class Policy> | |
114 | inline RealType mode(const non_central_f_distribution<RealType, Policy>& dist) | |
115 | { // mode. | |
116 | static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)"; | |
117 | ||
118 | RealType n = dist.degrees_of_freedom1(); | |
119 | RealType m = dist.degrees_of_freedom2(); | |
120 | RealType l = dist.non_centrality(); | |
121 | RealType r; | |
122 | if(!detail::check_df( | |
123 | function, | |
124 | n, &r, Policy()) | |
125 | || | |
126 | !detail::check_df( | |
127 | function, | |
128 | m, &r, Policy()) | |
129 | || | |
130 | !detail::check_non_centrality( | |
131 | function, | |
132 | l, | |
133 | &r, | |
134 | Policy())) | |
135 | return r; | |
136 | RealType guess = m > 2 ? RealType(m * (n + l) / (n * (m - 2))) : RealType(1); | |
137 | return detail::generic_find_mode( | |
138 | dist, | |
139 | guess, | |
140 | function); | |
141 | } | |
142 | ||
143 | template <class RealType, class Policy> | |
144 | inline RealType variance(const non_central_f_distribution<RealType, Policy>& dist) | |
145 | { // variance. | |
146 | const char* function = "variance(non_central_f_distribution<%1%> const&)"; | |
147 | RealType n = dist.degrees_of_freedom1(); | |
148 | RealType m = dist.degrees_of_freedom2(); | |
149 | RealType l = dist.non_centrality(); | |
150 | RealType r; | |
151 | if(!detail::check_df( | |
152 | function, | |
153 | n, &r, Policy()) | |
154 | || | |
155 | !detail::check_df( | |
156 | function, | |
157 | m, &r, Policy()) | |
158 | || | |
159 | !detail::check_non_centrality( | |
160 | function, | |
161 | l, | |
162 | &r, | |
163 | Policy())) | |
164 | return r; | |
165 | if(m <= 4) | |
166 | return policies::raise_domain_error( | |
167 | function, | |
168 | "Second degrees of freedom parameter was %1%, but must be > 4 !", | |
169 | m, Policy()); | |
170 | RealType result = 2 * m * m * ((n + l) * (n + l) | |
171 | + (m - 2) * (n + 2 * l)); | |
172 | result /= (m - 4) * (m - 2) * (m - 2) * n * n; | |
173 | return result; | |
174 | } | |
175 | ||
176 | // RealType standard_deviation(const non_central_f_distribution<RealType, Policy>& dist) | |
177 | // standard_deviation provided by derived accessors. | |
178 | ||
179 | template <class RealType, class Policy> | |
180 | inline RealType skewness(const non_central_f_distribution<RealType, Policy>& dist) | |
181 | { // skewness = sqrt(l). | |
182 | const char* function = "skewness(non_central_f_distribution<%1%> const&)"; | |
183 | BOOST_MATH_STD_USING | |
184 | RealType n = dist.degrees_of_freedom1(); | |
185 | RealType m = dist.degrees_of_freedom2(); | |
186 | RealType l = dist.non_centrality(); | |
187 | RealType r; | |
188 | if(!detail::check_df( | |
189 | function, | |
190 | n, &r, Policy()) | |
191 | || | |
192 | !detail::check_df( | |
193 | function, | |
194 | m, &r, Policy()) | |
195 | || | |
196 | !detail::check_non_centrality( | |
197 | function, | |
198 | l, | |
199 | &r, | |
200 | Policy())) | |
201 | return r; | |
202 | if(m <= 6) | |
203 | return policies::raise_domain_error( | |
204 | function, | |
205 | "Second degrees of freedom parameter was %1%, but must be > 6 !", | |
206 | m, Policy()); | |
207 | RealType result = 2 * constants::root_two<RealType>(); | |
208 | result *= sqrt(m - 4); | |
209 | result *= (n * (m + n - 2) *(m + 2 * n - 2) | |
210 | + 3 * (m + n - 2) * (m + 2 * n - 2) * l | |
211 | + 6 * (m + n - 2) * l * l + 2 * l * l * l); | |
212 | result /= (m - 6) * pow(n * (m + n - 2) + 2 * (m + n - 2) * l + l * l, RealType(1.5f)); | |
213 | return result; | |
214 | } | |
215 | ||
216 | template <class RealType, class Policy> | |
217 | inline RealType kurtosis_excess(const non_central_f_distribution<RealType, Policy>& dist) | |
218 | { | |
219 | const char* function = "kurtosis_excess(non_central_f_distribution<%1%> const&)"; | |
220 | BOOST_MATH_STD_USING | |
221 | RealType n = dist.degrees_of_freedom1(); | |
222 | RealType m = dist.degrees_of_freedom2(); | |
223 | RealType l = dist.non_centrality(); | |
224 | RealType r; | |
225 | if(!detail::check_df( | |
226 | function, | |
227 | n, &r, Policy()) | |
228 | || | |
229 | !detail::check_df( | |
230 | function, | |
231 | m, &r, Policy()) | |
232 | || | |
233 | !detail::check_non_centrality( | |
234 | function, | |
235 | l, | |
236 | &r, | |
237 | Policy())) | |
238 | return r; | |
239 | if(m <= 8) | |
240 | return policies::raise_domain_error( | |
241 | function, | |
242 | "Second degrees of freedom parameter was %1%, but must be > 8 !", | |
243 | m, Policy()); | |
244 | RealType l2 = l * l; | |
245 | RealType l3 = l2 * l; | |
246 | RealType l4 = l2 * l2; | |
247 | RealType result = (3 * (m - 4) * (n * (m + n - 2) | |
248 | * (4 * (m - 2) * (m - 2) | |
249 | + (m - 2) * (m + 10) * n | |
250 | + (10 + m) * n * n) | |
251 | + 4 * (m + n - 2) * (4 * (m - 2) * (m - 2) | |
252 | + (m - 2) * (10 + m) * n | |
253 | + (10 + m) * n * n) * l + 2 * (10 + m) | |
254 | * (m + n - 2) * (2 * m + 3 * n - 4) * l2 | |
255 | + 4 * (10 + m) * (-2 + m + n) * l3 | |
256 | + (10 + m) * l4)) | |
257 | / | |
258 | ((-8 + m) * (-6 + m) * boost::math::pow<2>(n * (-2 + m + n) | |
259 | + 2 * (-2 + m + n) * l + l2)); | |
260 | return result; | |
261 | } // kurtosis_excess | |
262 | ||
263 | template <class RealType, class Policy> | |
264 | inline RealType kurtosis(const non_central_f_distribution<RealType, Policy>& dist) | |
265 | { | |
266 | return kurtosis_excess(dist) + 3; | |
267 | } | |
268 | ||
269 | template <class RealType, class Policy> | |
270 | inline RealType pdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x) | |
271 | { // Probability Density/Mass Function. | |
272 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
273 | typedef typename policies::normalise< | |
274 | Policy, | |
275 | policies::promote_float<false>, | |
276 | policies::promote_double<false>, | |
277 | policies::discrete_quantile<>, | |
278 | policies::assert_undefined<> >::type forwarding_policy; | |
279 | ||
280 | value_type alpha = dist.degrees_of_freedom1() / 2; | |
281 | value_type beta = dist.degrees_of_freedom2() / 2; | |
282 | value_type y = x * alpha / beta; | |
283 | value_type r = pdf(boost::math::non_central_beta_distribution<value_type, forwarding_policy>(alpha, beta, dist.non_centrality()), y / (1 + y)); | |
284 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
285 | r * (dist.degrees_of_freedom1() / dist.degrees_of_freedom2()) / ((1 + y) * (1 + y)), | |
286 | "pdf(non_central_f_distribution<%1%>, %1%)"); | |
287 | ||
288 | ||
289 | template <class RealType, class Policy> | |
290 | RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x) | |
291 | { | |
292 | const char* function = "cdf(const non_central_f_distribution<%1%>&, %1%)"; | |
293 | RealType r; | |
294 | if(!detail::check_df( | |
295 | function, | |
296 | dist.degrees_of_freedom1(), &r, Policy()) | |
297 | || | |
298 | !detail::check_df( | |
299 | function, | |
300 | dist.degrees_of_freedom2(), &r, Policy()) | |
301 | || | |
302 | !detail::check_non_centrality( | |
303 | function, | |
304 | dist.non_centrality(), | |
305 | &r, | |
306 | Policy())) | |
307 | return r; | |
308 | ||
309 | if((x < 0) || !(boost::math::isfinite)(x)) | |
310 | { | |
311 | return policies::raise_domain_error<RealType>( | |
312 | function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy()); | |
313 | } | |
314 | ||
315 | RealType alpha = dist.degrees_of_freedom1() / 2; | |
316 | RealType beta = dist.degrees_of_freedom2() / 2; | |
317 | RealType y = x * alpha / beta; | |
318 | RealType c = y / (1 + y); | |
319 | RealType cp = 1 / (1 + y); | |
320 | // | |
321 | // To ensure accuracy, we pass both x and 1-x to the | |
322 | // non-central beta cdf routine, this ensures accuracy | |
323 | // even when we compute x to be ~ 1: | |
324 | // | |
325 | r = detail::non_central_beta_cdf(c, cp, alpha, beta, | |
326 | dist.non_centrality(), false, Policy()); | |
327 | return r; | |
328 | } // cdf | |
329 | ||
330 | template <class RealType, class Policy> | |
331 | RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c) | |
332 | { // Complemented Cumulative Distribution Function | |
333 | const char* function = "cdf(complement(const non_central_f_distribution<%1%>&, %1%))"; | |
334 | RealType r; | |
335 | if(!detail::check_df( | |
336 | function, | |
337 | c.dist.degrees_of_freedom1(), &r, Policy()) | |
338 | || | |
339 | !detail::check_df( | |
340 | function, | |
341 | c.dist.degrees_of_freedom2(), &r, Policy()) | |
342 | || | |
343 | !detail::check_non_centrality( | |
344 | function, | |
345 | c.dist.non_centrality(), | |
346 | &r, | |
347 | Policy())) | |
348 | return r; | |
349 | ||
350 | if((c.param < 0) || !(boost::math::isfinite)(c.param)) | |
351 | { | |
352 | return policies::raise_domain_error<RealType>( | |
353 | function, "Random Variable parameter was %1%, but must be > 0 !", c.param, Policy()); | |
354 | } | |
355 | ||
356 | RealType alpha = c.dist.degrees_of_freedom1() / 2; | |
357 | RealType beta = c.dist.degrees_of_freedom2() / 2; | |
358 | RealType y = c.param * alpha / beta; | |
359 | RealType x = y / (1 + y); | |
360 | RealType cx = 1 / (1 + y); | |
361 | // | |
362 | // To ensure accuracy, we pass both x and 1-x to the | |
363 | // non-central beta cdf routine, this ensures accuracy | |
364 | // even when we compute x to be ~ 1: | |
365 | // | |
366 | r = detail::non_central_beta_cdf(x, cx, alpha, beta, | |
367 | c.dist.non_centrality(), true, Policy()); | |
368 | return r; | |
369 | } // ccdf | |
370 | ||
371 | template <class RealType, class Policy> | |
372 | inline RealType quantile(const non_central_f_distribution<RealType, Policy>& dist, const RealType& p) | |
373 | { // Quantile (or Percent Point) function. | |
374 | RealType alpha = dist.degrees_of_freedom1() / 2; | |
375 | RealType beta = dist.degrees_of_freedom2() / 2; | |
376 | RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p); | |
377 | if(x == 1) | |
378 | return policies::raise_overflow_error<RealType>( | |
379 | "quantile(const non_central_f_distribution<%1%>&, %1%)", | |
380 | "Result of non central F quantile is too large to represent.", | |
381 | Policy()); | |
382 | return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1()); | |
383 | } // quantile | |
384 | ||
385 | template <class RealType, class Policy> | |
386 | inline RealType quantile(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c) | |
387 | { // Quantile (or Percent Point) function. | |
388 | RealType alpha = c.dist.degrees_of_freedom1() / 2; | |
389 | RealType beta = c.dist.degrees_of_freedom2() / 2; | |
390 | RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param)); | |
391 | if(x == 1) | |
392 | return policies::raise_overflow_error<RealType>( | |
393 | "quantile(complement(const non_central_f_distribution<%1%>&, %1%))", | |
394 | "Result of non central F quantile is too large to represent.", | |
395 | Policy()); | |
396 | return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1()); | |
397 | } // quantile complement. | |
398 | ||
399 | } // namespace math | |
400 | } // namespace boost | |
401 | ||
402 | // This include must be at the end, *after* the accessors | |
403 | // for this distribution have been defined, in order to | |
404 | // keep compilers that support two-phase lookup happy. | |
405 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
406 | ||
407 | #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP | |
408 | ||
409 | ||
410 |