]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // boost\math\distributions\non_central_beta.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_BETA_HPP | |
11 | #define BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP | |
12 | ||
13 | #include <boost/math/distributions/fwd.hpp> | |
14 | #include <boost/math/special_functions/beta.hpp> // for incomplete gamma. gamma_q | |
15 | #include <boost/math/distributions/complement.hpp> // complements | |
16 | #include <boost/math/distributions/beta.hpp> // central distribution | |
17 | #include <boost/math/distributions/detail/generic_mode.hpp> | |
18 | #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks | |
19 | #include <boost/math/special_functions/fpclassify.hpp> // isnan. | |
20 | #include <boost/math/tools/roots.hpp> // for root finding. | |
21 | #include <boost/math/tools/series.hpp> | |
22 | ||
23 | namespace boost | |
24 | { | |
25 | namespace math | |
26 | { | |
27 | ||
28 | template <class RealType, class Policy> | |
29 | class non_central_beta_distribution; | |
30 | ||
31 | namespace detail{ | |
32 | ||
33 | template <class T, class Policy> | |
34 | T non_central_beta_p(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0) | |
35 | { | |
36 | BOOST_MATH_STD_USING | |
37 | using namespace boost::math; | |
38 | // | |
39 | // Variables come first: | |
40 | // | |
41 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
42 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
43 | T l2 = lam / 2; | |
44 | // | |
45 | // k is the starting point for iteration, and is the | |
46 | // maximum of the poisson weighting term, | |
47 | // note that unlike other similar code, we do not set | |
48 | // k to zero, when l2 is small, as forward iteration | |
49 | // is unstable: | |
50 | // | |
51 | int k = itrunc(l2); | |
52 | if(k == 0) | |
53 | k = 1; | |
54 | // Starting Poisson weight: | |
55 | T pois = gamma_p_derivative(T(k+1), l2, pol); | |
56 | if(pois == 0) | |
57 | return init_val; | |
58 | // recurance term: | |
59 | T xterm; | |
60 | // Starting beta term: | |
61 | T beta = x < y | |
62 | ? detail::ibeta_imp(T(a + k), b, x, pol, false, true, &xterm) | |
63 | : detail::ibeta_imp(b, T(a + k), y, pol, true, true, &xterm); | |
64 | ||
65 | xterm *= y / (a + b + k - 1); | |
66 | T poisf(pois), betaf(beta), xtermf(xterm); | |
67 | T sum = init_val; | |
68 | ||
69 | if((beta == 0) && (xterm == 0)) | |
70 | return init_val; | |
71 | ||
72 | // | |
73 | // Backwards recursion first, this is the stable | |
74 | // direction for recursion: | |
75 | // | |
76 | T last_term = 0; | |
77 | boost::uintmax_t count = k; | |
78 | for(int i = k; i >= 0; --i) | |
79 | { | |
80 | T term = beta * pois; | |
81 | sum += term; | |
82 | if(((fabs(term/sum) < errtol) && (last_term >= term)) || (term == 0)) | |
83 | { | |
84 | count = k - i; | |
85 | break; | |
86 | } | |
87 | pois *= i / l2; | |
88 | beta += xterm; | |
89 | xterm *= (a + i - 1) / (x * (a + b + i - 2)); | |
90 | last_term = term; | |
91 | } | |
92 | for(int i = k + 1; ; ++i) | |
93 | { | |
94 | poisf *= l2 / i; | |
95 | xtermf *= (x * (a + b + i - 2)) / (a + i - 1); | |
96 | betaf -= xtermf; | |
97 | ||
98 | T term = poisf * betaf; | |
99 | sum += term; | |
100 | if((fabs(term/sum) < errtol) || (term == 0)) | |
101 | { | |
102 | break; | |
103 | } | |
104 | if(static_cast<boost::uintmax_t>(count + i - k) > max_iter) | |
105 | { | |
106 | return policies::raise_evaluation_error( | |
107 | "cdf(non_central_beta_distribution<%1%>, %1%)", | |
108 | "Series did not converge, closest value was %1%", sum, pol); | |
109 | } | |
110 | } | |
111 | return sum; | |
112 | } | |
113 | ||
114 | template <class T, class Policy> | |
115 | T non_central_beta_q(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0) | |
116 | { | |
117 | BOOST_MATH_STD_USING | |
118 | using namespace boost::math; | |
119 | // | |
120 | // Variables come first: | |
121 | // | |
122 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
123 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
124 | T l2 = lam / 2; | |
125 | // | |
126 | // k is the starting point for iteration, and is the | |
127 | // maximum of the poisson weighting term: | |
128 | // | |
129 | int k = itrunc(l2); | |
130 | T pois; | |
131 | if(k <= 30) | |
132 | { | |
133 | // | |
134 | // Might as well start at 0 since we'll likely have this number of terms anyway: | |
135 | // | |
136 | if(a + b > 1) | |
137 | k = 0; | |
138 | else if(k == 0) | |
139 | k = 1; | |
140 | } | |
141 | if(k == 0) | |
142 | { | |
143 | // Starting Poisson weight: | |
144 | pois = exp(-l2); | |
145 | } | |
146 | else | |
147 | { | |
148 | // Starting Poisson weight: | |
149 | pois = gamma_p_derivative(T(k+1), l2, pol); | |
150 | } | |
151 | if(pois == 0) | |
152 | return init_val; | |
153 | // recurance term: | |
154 | T xterm; | |
155 | // Starting beta term: | |
156 | T beta = x < y | |
157 | ? detail::ibeta_imp(T(a + k), b, x, pol, true, true, &xterm) | |
158 | : detail::ibeta_imp(b, T(a + k), y, pol, false, true, &xterm); | |
159 | ||
160 | xterm *= y / (a + b + k - 1); | |
161 | T poisf(pois), betaf(beta), xtermf(xterm); | |
162 | T sum = init_val; | |
163 | if((beta == 0) && (xterm == 0)) | |
164 | return init_val; | |
165 | // | |
166 | // Forwards recursion first, this is the stable | |
167 | // direction for recursion, and the location | |
168 | // of the bulk of the sum: | |
169 | // | |
170 | T last_term = 0; | |
171 | boost::uintmax_t count = 0; | |
172 | for(int i = k + 1; ; ++i) | |
173 | { | |
174 | poisf *= l2 / i; | |
175 | xtermf *= (x * (a + b + i - 2)) / (a + i - 1); | |
176 | betaf += xtermf; | |
177 | ||
178 | T term = poisf * betaf; | |
179 | sum += term; | |
180 | if((fabs(term/sum) < errtol) && (last_term >= term)) | |
181 | { | |
182 | count = i - k; | |
183 | break; | |
184 | } | |
185 | if(static_cast<boost::uintmax_t>(i - k) > max_iter) | |
186 | { | |
187 | return policies::raise_evaluation_error( | |
188 | "cdf(non_central_beta_distribution<%1%>, %1%)", | |
189 | "Series did not converge, closest value was %1%", sum, pol); | |
190 | } | |
191 | last_term = term; | |
192 | } | |
193 | for(int i = k; i >= 0; --i) | |
194 | { | |
195 | T term = beta * pois; | |
196 | sum += term; | |
197 | if(fabs(term/sum) < errtol) | |
198 | { | |
199 | break; | |
200 | } | |
201 | if(static_cast<boost::uintmax_t>(count + k - i) > max_iter) | |
202 | { | |
203 | return policies::raise_evaluation_error( | |
204 | "cdf(non_central_beta_distribution<%1%>, %1%)", | |
205 | "Series did not converge, closest value was %1%", sum, pol); | |
206 | } | |
207 | pois *= i / l2; | |
208 | beta -= xterm; | |
209 | xterm *= (a + i - 1) / (x * (a + b + i - 2)); | |
210 | } | |
211 | return sum; | |
212 | } | |
213 | ||
214 | template <class RealType, class Policy> | |
215 | inline RealType non_central_beta_cdf(RealType x, RealType y, RealType a, RealType b, RealType l, bool invert, const Policy&) | |
216 | { | |
217 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
218 | typedef typename policies::normalise< | |
219 | Policy, | |
220 | policies::promote_float<false>, | |
221 | policies::promote_double<false>, | |
222 | policies::discrete_quantile<>, | |
223 | policies::assert_undefined<> >::type forwarding_policy; | |
224 | ||
225 | BOOST_MATH_STD_USING | |
226 | ||
227 | if(x == 0) | |
228 | return invert ? 1.0f : 0.0f; | |
229 | if(y == 0) | |
230 | return invert ? 0.0f : 1.0f; | |
231 | value_type result; | |
232 | value_type c = a + b + l / 2; | |
233 | value_type cross = 1 - (b / c) * (1 + l / (2 * c * c)); | |
234 | if(l == 0) | |
235 | result = cdf(boost::math::beta_distribution<RealType, Policy>(a, b), x); | |
236 | else if(x > cross) | |
237 | { | |
238 | // Complement is the smaller of the two: | |
239 | result = detail::non_central_beta_q( | |
240 | static_cast<value_type>(a), | |
241 | static_cast<value_type>(b), | |
242 | static_cast<value_type>(l), | |
243 | static_cast<value_type>(x), | |
244 | static_cast<value_type>(y), | |
245 | forwarding_policy(), | |
246 | static_cast<value_type>(invert ? 0 : -1)); | |
247 | invert = !invert; | |
248 | } | |
249 | else | |
250 | { | |
251 | result = detail::non_central_beta_p( | |
252 | static_cast<value_type>(a), | |
253 | static_cast<value_type>(b), | |
254 | static_cast<value_type>(l), | |
255 | static_cast<value_type>(x), | |
256 | static_cast<value_type>(y), | |
257 | forwarding_policy(), | |
258 | static_cast<value_type>(invert ? -1 : 0)); | |
259 | } | |
260 | if(invert) | |
261 | result = -result; | |
262 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
263 | result, | |
264 | "boost::math::non_central_beta_cdf<%1%>(%1%, %1%, %1%)"); | |
265 | } | |
266 | ||
267 | template <class T, class Policy> | |
268 | struct nc_beta_quantile_functor | |
269 | { | |
270 | nc_beta_quantile_functor(const non_central_beta_distribution<T,Policy>& d, T t, bool c) | |
271 | : dist(d), target(t), comp(c) {} | |
272 | ||
273 | T operator()(const T& x) | |
274 | { | |
275 | return comp ? | |
276 | T(target - cdf(complement(dist, x))) | |
277 | : T(cdf(dist, x) - target); | |
278 | } | |
279 | ||
280 | private: | |
281 | non_central_beta_distribution<T,Policy> dist; | |
282 | T target; | |
283 | bool comp; | |
284 | }; | |
285 | ||
286 | // | |
287 | // This is more or less a copy of bracket_and_solve_root, but | |
288 | // modified to search only the interval [0,1] using similar | |
289 | // heuristics. | |
290 | // | |
291 | template <class F, class T, class Tol, class Policy> | |
292 | std::pair<T, T> bracket_and_solve_root_01(F f, const T& guess, T factor, bool rising, Tol tol, boost::uintmax_t& max_iter, const Policy& pol) | |
293 | { | |
294 | BOOST_MATH_STD_USING | |
295 | static const char* function = "boost::math::tools::bracket_and_solve_root_01<%1%>"; | |
296 | // | |
297 | // Set up inital brackets: | |
298 | // | |
299 | T a = guess; | |
300 | T b = a; | |
301 | T fa = f(a); | |
302 | T fb = fa; | |
303 | // | |
304 | // Set up invocation count: | |
305 | // | |
306 | boost::uintmax_t count = max_iter - 1; | |
307 | ||
308 | if((fa < 0) == (guess < 0 ? !rising : rising)) | |
309 | { | |
310 | // | |
311 | // Zero is to the right of b, so walk upwards | |
312 | // until we find it: | |
313 | // | |
314 | while((boost::math::sign)(fb) == (boost::math::sign)(fa)) | |
315 | { | |
316 | if(count == 0) | |
317 | { | |
318 | b = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol); | |
319 | return std::make_pair(a, b); | |
320 | } | |
321 | // | |
322 | // Heuristic: every 20 iterations we double the growth factor in case the | |
323 | // initial guess was *really* bad ! | |
324 | // | |
325 | if((max_iter - count) % 20 == 0) | |
326 | factor *= 2; | |
327 | // | |
328 | // Now go ahead and move are guess by "factor", | |
329 | // we do this by reducing 1-guess by factor: | |
330 | // | |
331 | a = b; | |
332 | fa = fb; | |
333 | b = 1 - ((1 - b) / factor); | |
334 | fb = f(b); | |
335 | --count; | |
336 | BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); | |
337 | } | |
338 | } | |
339 | else | |
340 | { | |
341 | // | |
342 | // Zero is to the left of a, so walk downwards | |
343 | // until we find it: | |
344 | // | |
345 | while((boost::math::sign)(fb) == (boost::math::sign)(fa)) | |
346 | { | |
347 | if(fabs(a) < tools::min_value<T>()) | |
348 | { | |
349 | // Escape route just in case the answer is zero! | |
350 | max_iter -= count; | |
351 | max_iter += 1; | |
352 | return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0)); | |
353 | } | |
354 | if(count == 0) | |
355 | { | |
356 | a = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol); | |
357 | return std::make_pair(a, b); | |
358 | } | |
359 | // | |
360 | // Heuristic: every 20 iterations we double the growth factor in case the | |
361 | // initial guess was *really* bad ! | |
362 | // | |
363 | if((max_iter - count) % 20 == 0) | |
364 | factor *= 2; | |
365 | // | |
366 | // Now go ahead and move are guess by "factor": | |
367 | // | |
368 | b = a; | |
369 | fb = fa; | |
370 | a /= factor; | |
371 | fa = f(a); | |
372 | --count; | |
373 | BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); | |
374 | } | |
375 | } | |
376 | max_iter -= count; | |
377 | max_iter += 1; | |
378 | std::pair<T, T> r = toms748_solve( | |
379 | f, | |
380 | (a < 0 ? b : a), | |
381 | (a < 0 ? a : b), | |
382 | (a < 0 ? fb : fa), | |
383 | (a < 0 ? fa : fb), | |
384 | tol, | |
385 | count, | |
386 | pol); | |
387 | max_iter += count; | |
388 | BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count); | |
389 | return r; | |
390 | } | |
391 | ||
392 | template <class RealType, class Policy> | |
393 | RealType nc_beta_quantile(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& p, bool comp) | |
394 | { | |
395 | static const char* function = "quantile(non_central_beta_distribution<%1%>, %1%)"; | |
396 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
397 | typedef typename policies::normalise< | |
398 | Policy, | |
399 | policies::promote_float<false>, | |
400 | policies::promote_double<false>, | |
401 | policies::discrete_quantile<>, | |
402 | policies::assert_undefined<> >::type forwarding_policy; | |
403 | ||
404 | value_type a = dist.alpha(); | |
405 | value_type b = dist.beta(); | |
406 | value_type l = dist.non_centrality(); | |
407 | value_type r; | |
408 | if(!beta_detail::check_alpha( | |
409 | function, | |
410 | a, &r, Policy()) | |
411 | || | |
412 | !beta_detail::check_beta( | |
413 | function, | |
414 | b, &r, Policy()) | |
415 | || | |
416 | !detail::check_non_centrality( | |
417 | function, | |
418 | l, | |
419 | &r, | |
420 | Policy()) | |
421 | || | |
422 | !detail::check_probability( | |
423 | function, | |
424 | static_cast<value_type>(p), | |
425 | &r, | |
426 | Policy())) | |
427 | return (RealType)r; | |
428 | // | |
429 | // Special cases first: | |
430 | // | |
431 | if(p == 0) | |
432 | return comp | |
433 | ? 1.0f | |
434 | : 0.0f; | |
435 | if(p == 1) | |
436 | return !comp | |
437 | ? 1.0f | |
438 | : 0.0f; | |
439 | ||
440 | value_type c = a + b + l / 2; | |
441 | value_type mean = 1 - (b / c) * (1 + l / (2 * c * c)); | |
442 | /* | |
443 | // | |
444 | // Calculate a normal approximation to the quantile, | |
445 | // uses mean and variance approximations from: | |
446 | // Algorithm AS 310: | |
447 | // Computing the Non-Central Beta Distribution Function | |
448 | // R. Chattamvelli; R. Shanmugam | |
449 | // Applied Statistics, Vol. 46, No. 1. (1997), pp. 146-156. | |
450 | // | |
451 | // Unfortunately, when this is wrong it tends to be *very* | |
452 | // wrong, so it's disabled for now, even though it often | |
453 | // gets the initial guess quite close. Probably we could | |
454 | // do much better by factoring in the skewness if only | |
455 | // we could calculate it.... | |
456 | // | |
457 | value_type delta = l / 2; | |
458 | value_type delta2 = delta * delta; | |
459 | value_type delta3 = delta * delta2; | |
460 | value_type delta4 = delta2 * delta2; | |
461 | value_type G = c * (c + 1) + delta; | |
462 | value_type alpha = a + b; | |
463 | value_type alpha2 = alpha * alpha; | |
464 | value_type eta = (2 * alpha + 1) * (2 * alpha + 1) + 1; | |
465 | value_type H = 3 * alpha2 + 5 * alpha + 2; | |
466 | value_type F = alpha2 * (alpha + 1) + H * delta | |
467 | + (2 * alpha + 4) * delta2 + delta3; | |
468 | value_type P = (3 * alpha + 1) * (9 * alpha + 17) | |
469 | + 2 * alpha * (3 * alpha + 2) * (3 * alpha + 4) + 15; | |
470 | value_type Q = 54 * alpha2 + 162 * alpha + 130; | |
471 | value_type R = 6 * (6 * alpha + 11); | |
472 | value_type D = delta | |
473 | * (H * H + 2 * P * delta + Q * delta2 + R * delta3 + 9 * delta4); | |
474 | value_type variance = (b / G) | |
475 | * (1 + delta * (l * l + 3 * l + eta) / (G * G)) | |
476 | - (b * b / F) * (1 + D / (F * F)); | |
477 | value_type sd = sqrt(variance); | |
478 | ||
479 | value_type guess = comp | |
480 | ? quantile(complement(normal_distribution<RealType, Policy>(static_cast<RealType>(mean), static_cast<RealType>(sd)), p)) | |
481 | : quantile(normal_distribution<RealType, Policy>(static_cast<RealType>(mean), static_cast<RealType>(sd)), p); | |
482 | ||
483 | if(guess >= 1) | |
484 | guess = mean; | |
485 | if(guess <= tools::min_value<value_type>()) | |
486 | guess = mean; | |
487 | */ | |
488 | value_type guess = mean; | |
489 | detail::nc_beta_quantile_functor<value_type, Policy> | |
490 | f(non_central_beta_distribution<value_type, Policy>(a, b, l), p, comp); | |
491 | tools::eps_tolerance<value_type> tol(policies::digits<RealType, Policy>()); | |
492 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
493 | ||
494 | std::pair<value_type, value_type> ir | |
495 | = bracket_and_solve_root_01( | |
496 | f, guess, value_type(2.5), true, tol, | |
497 | max_iter, Policy()); | |
498 | value_type result = ir.first + (ir.second - ir.first) / 2; | |
499 | ||
500 | if(max_iter >= policies::get_max_root_iterations<Policy>()) | |
501 | { | |
502 | return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" | |
503 | " either there is no answer to quantile of the non central beta distribution" | |
504 | " or the answer is infinite. Current best guess is %1%", | |
505 | policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
506 | result, | |
507 | function), Policy()); | |
508 | } | |
509 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
510 | result, | |
511 | function); | |
512 | } | |
513 | ||
514 | template <class T, class Policy> | |
515 | T non_central_beta_pdf(T a, T b, T lam, T x, T y, const Policy& pol) | |
516 | { | |
517 | BOOST_MATH_STD_USING | |
518 | // | |
519 | // Special cases: | |
520 | // | |
521 | if((x == 0) || (y == 0)) | |
522 | return 0; | |
523 | // | |
524 | // Variables come first: | |
525 | // | |
526 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
527 | T errtol = boost::math::policies::get_epsilon<T, Policy>(); | |
528 | T l2 = lam / 2; | |
529 | // | |
530 | // k is the starting point for iteration, and is the | |
531 | // maximum of the poisson weighting term: | |
532 | // | |
533 | int k = itrunc(l2); | |
534 | // Starting Poisson weight: | |
535 | T pois = gamma_p_derivative(T(k+1), l2, pol); | |
536 | // Starting beta term: | |
537 | T beta = x < y ? | |
538 | ibeta_derivative(a + k, b, x, pol) | |
539 | : ibeta_derivative(b, a + k, y, pol); | |
540 | T sum = 0; | |
541 | T poisf(pois); | |
542 | T betaf(beta); | |
543 | ||
544 | // | |
545 | // Stable backwards recursion first: | |
546 | // | |
547 | boost::uintmax_t count = k; | |
548 | for(int i = k; i >= 0; --i) | |
549 | { | |
550 | T term = beta * pois; | |
551 | sum += term; | |
552 | if((fabs(term/sum) < errtol) || (term == 0)) | |
553 | { | |
554 | count = k - i; | |
555 | break; | |
556 | } | |
557 | pois *= i / l2; | |
558 | beta *= (a + i - 1) / (x * (a + i + b - 1)); | |
559 | } | |
560 | for(int i = k + 1; ; ++i) | |
561 | { | |
562 | poisf *= l2 / i; | |
563 | betaf *= x * (a + b + i - 1) / (a + i - 1); | |
564 | ||
565 | T term = poisf * betaf; | |
566 | sum += term; | |
567 | if((fabs(term/sum) < errtol) || (term == 0)) | |
568 | { | |
569 | break; | |
570 | } | |
571 | if(static_cast<boost::uintmax_t>(count + i - k) > max_iter) | |
572 | { | |
573 | return policies::raise_evaluation_error( | |
574 | "pdf(non_central_beta_distribution<%1%>, %1%)", | |
575 | "Series did not converge, closest value was %1%", sum, pol); | |
576 | } | |
577 | } | |
578 | return sum; | |
579 | } | |
580 | ||
581 | template <class RealType, class Policy> | |
582 | RealType nc_beta_pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x) | |
583 | { | |
584 | BOOST_MATH_STD_USING | |
585 | static const char* function = "pdf(non_central_beta_distribution<%1%>, %1%)"; | |
586 | typedef typename policies::evaluation<RealType, Policy>::type value_type; | |
587 | typedef typename policies::normalise< | |
588 | Policy, | |
589 | policies::promote_float<false>, | |
590 | policies::promote_double<false>, | |
591 | policies::discrete_quantile<>, | |
592 | policies::assert_undefined<> >::type forwarding_policy; | |
593 | ||
594 | value_type a = dist.alpha(); | |
595 | value_type b = dist.beta(); | |
596 | value_type l = dist.non_centrality(); | |
597 | value_type r; | |
598 | if(!beta_detail::check_alpha( | |
599 | function, | |
600 | a, &r, Policy()) | |
601 | || | |
602 | !beta_detail::check_beta( | |
603 | function, | |
604 | b, &r, Policy()) | |
605 | || | |
606 | !detail::check_non_centrality( | |
607 | function, | |
608 | l, | |
609 | &r, | |
610 | Policy()) | |
611 | || | |
612 | !beta_detail::check_x( | |
613 | function, | |
614 | static_cast<value_type>(x), | |
615 | &r, | |
616 | Policy())) | |
617 | return (RealType)r; | |
618 | ||
619 | if(l == 0) | |
620 | return pdf(boost::math::beta_distribution<RealType, Policy>(dist.alpha(), dist.beta()), x); | |
621 | return policies::checked_narrowing_cast<RealType, forwarding_policy>( | |
622 | non_central_beta_pdf(a, b, l, static_cast<value_type>(x), value_type(1 - static_cast<value_type>(x)), forwarding_policy()), | |
623 | "function"); | |
624 | } | |
625 | ||
626 | template <class T> | |
627 | struct hypergeometric_2F2_sum | |
628 | { | |
629 | typedef T result_type; | |
630 | hypergeometric_2F2_sum(T a1_, T a2_, T b1_, T b2_, T z_) : a1(a1_), a2(a2_), b1(b1_), b2(b2_), z(z_), term(1), k(0) {} | |
631 | T operator()() | |
632 | { | |
633 | T result = term; | |
634 | term *= a1 * a2 / (b1 * b2); | |
635 | a1 += 1; | |
636 | a2 += 1; | |
637 | b1 += 1; | |
638 | b2 += 1; | |
639 | k += 1; | |
640 | term /= k; | |
641 | term *= z; | |
642 | return result; | |
643 | } | |
644 | T a1, a2, b1, b2, z, term, k; | |
645 | }; | |
646 | ||
647 | template <class T, class Policy> | |
648 | T hypergeometric_2F2(T a1, T a2, T b1, T b2, T z, const Policy& pol) | |
649 | { | |
650 | typedef typename policies::evaluation<T, Policy>::type value_type; | |
651 | ||
652 | const char* function = "boost::math::detail::hypergeometric_2F2<%1%>(%1%,%1%,%1%,%1%,%1%)"; | |
653 | ||
654 | hypergeometric_2F2_sum<value_type> s(a1, a2, b1, b2, z); | |
655 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
656 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
657 | value_type zero = 0; | |
658 | value_type result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<value_type, Policy>(), max_iter, zero); | |
659 | #else | |
660 | value_type result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<value_type, Policy>(), max_iter); | |
661 | #endif | |
662 | policies::check_series_iterations<T>(function, max_iter, pol); | |
663 | return policies::checked_narrowing_cast<T, Policy>(result, function); | |
664 | } | |
665 | ||
666 | } // namespace detail | |
667 | ||
668 | template <class RealType = double, class Policy = policies::policy<> > | |
669 | class non_central_beta_distribution | |
670 | { | |
671 | public: | |
672 | typedef RealType value_type; | |
673 | typedef Policy policy_type; | |
674 | ||
675 | non_central_beta_distribution(RealType a_, RealType b_, RealType lambda) : a(a_), b(b_), ncp(lambda) | |
676 | { | |
677 | const char* function = "boost::math::non_central_beta_distribution<%1%>::non_central_beta_distribution(%1%,%1%)"; | |
678 | RealType r; | |
679 | beta_detail::check_alpha( | |
680 | function, | |
681 | a, &r, Policy()); | |
682 | beta_detail::check_beta( | |
683 | function, | |
684 | b, &r, Policy()); | |
685 | detail::check_non_centrality( | |
686 | function, | |
687 | lambda, | |
688 | &r, | |
689 | Policy()); | |
690 | } // non_central_beta_distribution constructor. | |
691 | ||
692 | RealType alpha() const | |
693 | { // Private data getter function. | |
694 | return a; | |
695 | } | |
696 | RealType beta() const | |
697 | { // Private data getter function. | |
698 | return b; | |
699 | } | |
700 | RealType non_centrality() const | |
701 | { // Private data getter function. | |
702 | return ncp; | |
703 | } | |
704 | private: | |
705 | // Data member, initialized by constructor. | |
706 | RealType a; // alpha. | |
707 | RealType b; // beta. | |
708 | RealType ncp; // non-centrality parameter | |
709 | }; // template <class RealType, class Policy> class non_central_beta_distribution | |
710 | ||
711 | typedef non_central_beta_distribution<double> non_central_beta; // Reserved name of type double. | |
712 | ||
713 | // Non-member functions to give properties of the distribution. | |
714 | ||
715 | template <class RealType, class Policy> | |
716 | inline const std::pair<RealType, RealType> range(const non_central_beta_distribution<RealType, Policy>& /* dist */) | |
717 | { // Range of permissible values for random variable k. | |
718 | using boost::math::tools::max_value; | |
719 | return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); | |
720 | } | |
721 | ||
722 | template <class RealType, class Policy> | |
723 | inline const std::pair<RealType, RealType> support(const non_central_beta_distribution<RealType, Policy>& /* dist */) | |
724 | { // Range of supported values for random variable k. | |
725 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
726 | using boost::math::tools::max_value; | |
727 | return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); | |
728 | } | |
729 | ||
730 | template <class RealType, class Policy> | |
731 | inline RealType mode(const non_central_beta_distribution<RealType, Policy>& dist) | |
732 | { // mode. | |
733 | static const char* function = "mode(non_central_beta_distribution<%1%> const&)"; | |
734 | ||
735 | RealType a = dist.alpha(); | |
736 | RealType b = dist.beta(); | |
737 | RealType l = dist.non_centrality(); | |
738 | RealType r; | |
739 | if(!beta_detail::check_alpha( | |
740 | function, | |
741 | a, &r, Policy()) | |
742 | || | |
743 | !beta_detail::check_beta( | |
744 | function, | |
745 | b, &r, Policy()) | |
746 | || | |
747 | !detail::check_non_centrality( | |
748 | function, | |
749 | l, | |
750 | &r, | |
751 | Policy())) | |
752 | return (RealType)r; | |
753 | RealType c = a + b + l / 2; | |
754 | RealType mean = 1 - (b / c) * (1 + l / (2 * c * c)); | |
755 | return detail::generic_find_mode_01( | |
756 | dist, | |
757 | mean, | |
758 | function); | |
759 | } | |
760 | ||
761 | // | |
762 | // We don't have the necessary information to implement | |
763 | // these at present. These are just disabled for now, | |
764 | // prototypes retained so we can fill in the blanks | |
765 | // later: | |
766 | // | |
767 | template <class RealType, class Policy> | |
768 | inline RealType mean(const non_central_beta_distribution<RealType, Policy>& dist) | |
769 | { | |
770 | BOOST_MATH_STD_USING | |
771 | RealType a = dist.alpha(); | |
772 | RealType b = dist.beta(); | |
773 | RealType d = dist.non_centrality(); | |
774 | RealType apb = a + b; | |
775 | return exp(-d / 2) * a * detail::hypergeometric_2F2<RealType, Policy>(1 + a, apb, a, 1 + apb, d / 2, Policy()) / apb; | |
776 | } // mean | |
777 | ||
778 | template <class RealType, class Policy> | |
779 | inline RealType variance(const non_central_beta_distribution<RealType, Policy>& dist) | |
780 | { | |
781 | // | |
782 | // Relative error of this function may be arbitarily large... absolute | |
783 | // error will be small however... that's the best we can do for now. | |
784 | // | |
785 | BOOST_MATH_STD_USING | |
786 | RealType a = dist.alpha(); | |
787 | RealType b = dist.beta(); | |
788 | RealType d = dist.non_centrality(); | |
789 | RealType apb = a + b; | |
790 | RealType result = detail::hypergeometric_2F2(RealType(1 + a), apb, a, RealType(1 + apb), RealType(d / 2), Policy()); | |
791 | result *= result * -exp(-d) * a * a / (apb * apb); | |
792 | result += exp(-d / 2) * a * (1 + a) * detail::hypergeometric_2F2(RealType(2 + a), apb, a, RealType(2 + apb), RealType(d / 2), Policy()) / (apb * (1 + apb)); | |
793 | return result; | |
794 | } | |
795 | ||
796 | // RealType standard_deviation(const non_central_beta_distribution<RealType, Policy>& dist) | |
797 | // standard_deviation provided by derived accessors. | |
798 | template <class RealType, class Policy> | |
799 | inline RealType skewness(const non_central_beta_distribution<RealType, Policy>& /*dist*/) | |
800 | { // skewness = sqrt(l). | |
801 | const char* function = "boost::math::non_central_beta_distribution<%1%>::skewness()"; | |
802 | typedef typename Policy::assert_undefined_type assert_type; | |
803 | BOOST_STATIC_ASSERT(assert_type::value == 0); | |
804 | ||
805 | return policies::raise_evaluation_error<RealType>( | |
806 | function, | |
807 | "This function is not yet implemented, the only sensible result is %1%.", | |
808 | std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? | |
809 | } | |
810 | ||
811 | template <class RealType, class Policy> | |
812 | inline RealType kurtosis_excess(const non_central_beta_distribution<RealType, Policy>& /*dist*/) | |
813 | { | |
814 | const char* function = "boost::math::non_central_beta_distribution<%1%>::kurtosis_excess()"; | |
815 | typedef typename Policy::assert_undefined_type assert_type; | |
816 | BOOST_STATIC_ASSERT(assert_type::value == 0); | |
817 | ||
818 | return policies::raise_evaluation_error<RealType>( | |
819 | function, | |
820 | "This function is not yet implemented, the only sensible result is %1%.", | |
821 | std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? | |
822 | } // kurtosis_excess | |
823 | ||
824 | template <class RealType, class Policy> | |
825 | inline RealType kurtosis(const non_central_beta_distribution<RealType, Policy>& dist) | |
826 | { | |
827 | return kurtosis_excess(dist) + 3; | |
828 | } | |
829 | ||
830 | template <class RealType, class Policy> | |
831 | inline RealType pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x) | |
832 | { // Probability Density/Mass Function. | |
833 | return detail::nc_beta_pdf(dist, x); | |
834 | ||
835 | ||
836 | template <class RealType, class Policy> | |
837 | RealType cdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x) | |
838 | { | |
839 | const char* function = "boost::math::non_central_beta_distribution<%1%>::cdf(%1%)"; | |
840 | RealType a = dist.alpha(); | |
841 | RealType b = dist.beta(); | |
842 | RealType l = dist.non_centrality(); | |
843 | RealType r; | |
844 | if(!beta_detail::check_alpha( | |
845 | function, | |
846 | a, &r, Policy()) | |
847 | || | |
848 | !beta_detail::check_beta( | |
849 | function, | |
850 | b, &r, Policy()) | |
851 | || | |
852 | !detail::check_non_centrality( | |
853 | function, | |
854 | l, | |
855 | &r, | |
856 | Policy()) | |
857 | || | |
858 | !beta_detail::check_x( | |
859 | function, | |
860 | x, | |
861 | &r, | |
862 | Policy())) | |
863 | return (RealType)r; | |
864 | ||
865 | if(l == 0) | |
866 | return cdf(beta_distribution<RealType, Policy>(a, b), x); | |
867 | ||
868 | return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, false, Policy()); | |
869 | } // cdf | |
870 | ||
871 | template <class RealType, class Policy> | |
872 | RealType cdf(const complemented2_type<non_central_beta_distribution<RealType, Policy>, RealType>& c) | |
873 | { // Complemented Cumulative Distribution Function | |
874 | const char* function = "boost::math::non_central_beta_distribution<%1%>::cdf(%1%)"; | |
875 | non_central_beta_distribution<RealType, Policy> const& dist = c.dist; | |
876 | RealType a = dist.alpha(); | |
877 | RealType b = dist.beta(); | |
878 | RealType l = dist.non_centrality(); | |
879 | RealType x = c.param; | |
880 | RealType r; | |
881 | if(!beta_detail::check_alpha( | |
882 | function, | |
883 | a, &r, Policy()) | |
884 | || | |
885 | !beta_detail::check_beta( | |
886 | function, | |
887 | b, &r, Policy()) | |
888 | || | |
889 | !detail::check_non_centrality( | |
890 | function, | |
891 | l, | |
892 | &r, | |
893 | Policy()) | |
894 | || | |
895 | !beta_detail::check_x( | |
896 | function, | |
897 | x, | |
898 | &r, | |
899 | Policy())) | |
900 | return (RealType)r; | |
901 | ||
902 | if(l == 0) | |
903 | return cdf(complement(beta_distribution<RealType, Policy>(a, b), x)); | |
904 | ||
905 | return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, true, Policy()); | |
906 | } // ccdf | |
907 | ||
908 | template <class RealType, class Policy> | |
909 | inline RealType quantile(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& p) | |
910 | { // Quantile (or Percent Point) function. | |
911 | return detail::nc_beta_quantile(dist, p, false); | |
912 | } // quantile | |
913 | ||
914 | template <class RealType, class Policy> | |
915 | inline RealType quantile(const complemented2_type<non_central_beta_distribution<RealType, Policy>, RealType>& c) | |
916 | { // Quantile (or Percent Point) function. | |
917 | return detail::nc_beta_quantile(c.dist, c.param, true); | |
918 | } // quantile complement. | |
919 | ||
920 | } // namespace math | |
921 | } // namespace boost | |
922 | ||
923 | // This include must be at the end, *after* the accessors | |
924 | // for this distribution have been defined, in order to | |
925 | // keep compilers that support two-phase lookup happy. | |
926 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
927 | ||
928 | #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP | |
929 |