1 // (C) Copyright John Maddock 2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 #ifndef BOOST_MATH_SPECIAL_BETA_HPP
7 #define BOOST_MATH_SPECIAL_BETA_HPP
13 #include <boost/math/special_functions/math_fwd.hpp>
14 #include <boost/math/tools/config.hpp>
15 #include <boost/math/special_functions/gamma.hpp>
16 #include <boost/math/special_functions/factorials.hpp>
17 #include <boost/math/special_functions/erf.hpp>
18 #include <boost/math/special_functions/log1p.hpp>
19 #include <boost/math/special_functions/expm1.hpp>
20 #include <boost/math/special_functions/trunc.hpp>
21 #include <boost/math/tools/roots.hpp>
22 #include <boost/static_assert.hpp>
23 #include <boost/config/no_tr1/cmath.hpp>
25 namespace boost{ namespace math{
30 // Implementation of Beta(a,b) using the Lanczos approximation:
32 template <class T, class Lanczos, class Policy>
33 T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
35 BOOST_MATH_STD_USING // for ADL of std names
38 policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol);
40 policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol);
48 if((c == a) && (b < tools::epsilon<T>()))
49 return boost::math::tgamma(b, pol);
50 else if((c == b) && (a < tools::epsilon<T>()))
51 return boost::math::tgamma(a, pol);
59 // This code appears to be no longer necessary: it was
60 // used to offset errors introduced from the Lanczos
61 // approximation, but the current Lanczos approximations
62 // are sufficiently accurate for all z that we can ditch
63 // this. It remains in the file for future reference...
65 // If a or b are less than 1, shift to greater than 1:
83 // Lanczos calculation:
84 T agh = a + Lanczos::g() - T(0.5);
85 T bgh = b + Lanczos::g() - T(0.5);
86 T cgh = c + Lanczos::g() - T(0.5);
87 result = Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c);
88 T ambh = a - T(0.5) - b;
89 if((fabs(b * ambh) < (cgh * 100)) && (a > 100))
91 // Special case where the base of the power term is close to 1
92 // compute (1+x)^y instead:
93 result *= exp(ambh * boost::math::log1p(-b / cgh, pol));
97 result *= pow(agh / cgh, a - T(0.5) - b);
100 // this avoids possible overflow, but appears to be marginally less accurate:
101 result *= pow((agh / cgh) * (bgh / cgh), b);
103 result *= pow((agh * bgh) / (cgh * cgh), b);
104 result *= sqrt(boost::math::constants::e<T>() / bgh);
106 // If a and b were originally less than 1 we need to scale the result:
110 } // template <class T, class Lanczos> beta_imp(T a, T b, const Lanczos&)
113 // Generic implementation of Beta(a,b) without Lanczos approximation support
114 // (Caution this is slow!!!):
116 template <class T, class Policy>
117 T beta_imp(T a, T b, const lanczos::undefined_lanczos& /* l */, const Policy& pol)
122 policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol);
124 policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol);
132 if((c == a) && (b < tools::epsilon<T>()))
133 return boost::math::tgamma(b, pol);
134 else if((c == b) && (a < tools::epsilon<T>()))
135 return boost::math::tgamma(a, pol);
141 // shift to a and b > 1 if required:
157 // set integration limits:
158 T la = (std::max)(T(10), a);
159 T lb = (std::max)(T(10), b);
160 T lc = (std::max)(T(10), T(a+b));
162 // calculate the fraction parts:
163 T sa = detail::lower_gamma_series(a, la, pol) / a;
164 sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
165 T sb = detail::lower_gamma_series(b, lb, pol) / b;
166 sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
167 T sc = detail::lower_gamma_series(c, lc, pol) / c;
168 sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
170 // and the exponent part:
171 result = exp(lc - la - lb) * pow(la/lc, a) * pow(lb/lc, b);
174 result *= sa * sb / sc;
176 // if a and b were originally less than 1 we need to scale the result:
180 } // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l)
184 // Compute the leading power terms in the incomplete Beta:
186 // (x^a)(y^b)/Beta(a,b) when normalised, and
187 // (x^a)(y^b) otherwise.
189 // Almost all of the error in the incomplete beta comes from this
190 // function: particularly when a and b are large. Computing large
191 // powers are *hard* though, and using logarithms just leads to
192 // horrendous cancellation errors.
194 template <class T, class Lanczos, class Policy>
195 T ibeta_power_terms(T a,
207 // can we do better here?
208 return pow(x, a) * pow(y, b);
216 // combine power terms with Lanczos approximation:
217 T agh = a + Lanczos::g() - T(0.5);
218 T bgh = b + Lanczos::g() - T(0.5);
219 T cgh = c + Lanczos::g() - T(0.5);
220 result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
222 // l1 and l2 are the base of the exponents minus one:
223 T l1 = (x * b - y * agh) / agh;
224 T l2 = (y * a - x * bgh) / bgh;
225 if(((std::min)(fabs(l1), fabs(l2)) < 0.2))
227 // when the base of the exponent is very near 1 we get really
228 // gross errors unless extra care is taken:
229 if((l1 * l2 > 0) || ((std::min)(a, b) < 1))
232 // This first branch handles the simple cases where either:
234 // * The two power terms both go in the same direction
235 // (towards zero or towards infinity). In this case if either
236 // term overflows or underflows, then the product of the two must
238 // *Alternatively if one exponent is less than one, then we
239 // can't productively use it to eliminate overflow or underflow
240 // from the other term. Problems with spurious overflow/underflow
241 // can't be ruled out in this case, but it is *very* unlikely
242 // since one of the power terms will evaluate to a number close to 1.
246 result *= exp(a * boost::math::log1p(l1, pol));
247 BOOST_MATH_INSTRUMENT_VARIABLE(result);
251 result *= pow((x * cgh) / agh, a);
252 BOOST_MATH_INSTRUMENT_VARIABLE(result);
256 result *= exp(b * boost::math::log1p(l2, pol));
257 BOOST_MATH_INSTRUMENT_VARIABLE(result);
261 result *= pow((y * cgh) / bgh, b);
262 BOOST_MATH_INSTRUMENT_VARIABLE(result);
265 else if((std::max)(fabs(l1), fabs(l2)) < 0.5)
268 // Both exponents are near one and both the exponents are
269 // greater than one and further these two
270 // power terms tend in opposite directions (one towards zero,
271 // the other towards infinity), so we have to combine the terms
272 // to avoid any risk of overflow or underflow.
274 // We do this by moving one power term inside the other, we have:
276 // (1 + l1)^a * (1 + l2)^b
277 // = ((1 + l1)*(1 + l2)^(b/a))^a
278 // = (1 + l1 + l3 + l1*l3)^a ; l3 = (1 + l2)^(b/a) - 1
279 // = exp((b/a) * log(1 + l2)) - 1
281 // The tricky bit is deciding which term to move inside :-)
282 // By preference we move the larger term inside, so that the
283 // size of the largest exponent is reduced. However, that can
284 // only be done as long as l3 (see above) is also small.
286 bool small_a = a < b;
288 if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1)))
290 T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol);
291 l3 = l1 + l3 + l3 * l1;
292 l3 = a * boost::math::log1p(l3, pol);
294 BOOST_MATH_INSTRUMENT_VARIABLE(result);
298 T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol);
299 l3 = l2 + l3 + l3 * l2;
300 l3 = b * boost::math::log1p(l3, pol);
302 BOOST_MATH_INSTRUMENT_VARIABLE(result);
305 else if(fabs(l1) < fabs(l2))
307 // First base near 1 only:
308 T l = a * boost::math::log1p(l1, pol)
309 + b * log((y * cgh) / bgh);
311 BOOST_MATH_INSTRUMENT_VARIABLE(result);
315 // Second base near 1 only:
316 T l = b * boost::math::log1p(l2, pol)
317 + a * log((x * cgh) / agh);
319 BOOST_MATH_INSTRUMENT_VARIABLE(result);
325 T b1 = (x * cgh) / agh;
326 T b2 = (y * cgh) / bgh;
329 BOOST_MATH_INSTRUMENT_VARIABLE(b1);
330 BOOST_MATH_INSTRUMENT_VARIABLE(b2);
331 BOOST_MATH_INSTRUMENT_VARIABLE(l1);
332 BOOST_MATH_INSTRUMENT_VARIABLE(l2);
333 if((l1 >= tools::log_max_value<T>())
334 || (l1 <= tools::log_min_value<T>())
335 || (l2 >= tools::log_max_value<T>())
336 || (l2 <= tools::log_min_value<T>())
339 // Oops, overflow, sidestep:
341 result *= pow(pow(b2, b/a) * b1, a);
343 result *= pow(pow(b1, a/b) * b2, b);
344 BOOST_MATH_INSTRUMENT_VARIABLE(result);
348 // finally the normal case:
349 result *= pow(b1, a) * pow(b2, b);
350 BOOST_MATH_INSTRUMENT_VARIABLE(result);
353 // combine with the leftover terms from the Lanczos approximation:
354 result *= sqrt(bgh / boost::math::constants::e<T>());
355 result *= sqrt(agh / cgh);
358 BOOST_MATH_INSTRUMENT_VARIABLE(result);
363 // Compute the leading power terms in the incomplete Beta:
365 // (x^a)(y^b)/Beta(a,b) when normalised, and
366 // (x^a)(y^b) otherwise.
368 // Almost all of the error in the incomplete beta comes from this
369 // function: particularly when a and b are large. Computing large
370 // powers are *hard* though, and using logarithms just leads to
371 // horrendous cancellation errors.
373 // This version is generic, slow, and does not use the Lanczos approximation.
375 template <class T, class Policy>
376 T ibeta_power_terms(T a,
380 const boost::math::lanczos::undefined_lanczos&,
388 return pow(x, a) * pow(y, b);
391 T result= 0; // assignment here silences warnings later
395 // integration limits for the gamma functions:
396 //T la = (std::max)(T(10), a);
397 //T lb = (std::max)(T(10), b);
398 //T lc = (std::max)(T(10), a+b);
402 // gamma function partials:
403 T sa = detail::lower_gamma_series(a, la, pol) / a;
404 sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
405 T sb = detail::lower_gamma_series(b, lb, pol) / b;
406 sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
407 T sc = detail::lower_gamma_series(c, lc, pol) / c;
408 sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
409 // gamma function powers combined with incomplete beta powers:
411 T b1 = (x * lc) / la;
412 T b2 = (y * lc) / lb;
417 if((lb1 >= tools::log_max_value<T>())
418 || (lb1 <= tools::log_min_value<T>())
419 || (lb2 >= tools::log_max_value<T>())
420 || (lb2 <= tools::log_min_value<T>())
421 || (e1 >= tools::log_max_value<T>())
422 || (e1 <= tools::log_min_value<T>())
425 result = exp(lb1 + lb2 - e1);
430 if((fabs(b1 - 1) * a < 10) && (a > 1))
431 p1 = exp(a * boost::math::log1p((x * b - y * la) / la, pol));
434 if((fabs(b2 - 1) * b < 10) && (b > 1))
435 p2 = exp(b * boost::math::log1p((y * a - x * lb) / lb, pol));
439 result = p1 * p2 / p3;
441 // and combine with the remaining gamma function components:
442 result /= sa * sb / sc;
447 // Series approximation to the incomplete beta:
450 struct ibeta_series_t
452 typedef T result_type;
453 ibeta_series_t(T a_, T b_, T x_, T mult) : result(mult), x(x_), apn(a_), poch(1-b_), n(1) {}
458 result *= poch * x / n;
464 T result, x, apn, poch;
468 template <class T, class Lanczos, class Policy>
469 T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol)
475 BOOST_ASSERT((p_derivative == 0) || normalised);
481 // incomplete beta power term, combined with the Lanczos approximation:
482 T agh = a + Lanczos::g() - T(0.5);
483 T bgh = b + Lanczos::g() - T(0.5);
484 T cgh = c + Lanczos::g() - T(0.5);
485 result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
487 result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol));
489 result *= pow(cgh / bgh, b - 0.5f);
490 result *= pow(x * cgh / agh, a);
491 result *= sqrt(agh / boost::math::constants::e<T>());
495 *p_derivative = result * pow(y, b);
496 BOOST_ASSERT(*p_derivative >= 0);
501 // Non-normalised, just compute the power:
504 if(result < tools::min_value<T>())
505 return s0; // Safeguard: series can't cope with denorms.
506 ibeta_series_t<T> s(a, b, x, result);
507 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
508 result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
509 policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol);
513 // Incomplete Beta series again, this time without Lanczos support:
515 template <class T, class Policy>
516 T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol)
521 BOOST_ASSERT((p_derivative == 0) || normalised);
527 // figure out integration limits for the gamma function:
528 //T la = (std::max)(T(10), a);
529 //T lb = (std::max)(T(10), b);
530 //T lc = (std::max)(T(10), a+b);
535 // calculate the gamma parts:
536 T sa = detail::lower_gamma_series(a, la, pol) / a;
537 sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
538 T sb = detail::lower_gamma_series(b, lb, pol) / b;
539 sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
540 T sc = detail::lower_gamma_series(c, lc, pol) / c;
541 sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
543 // and their combined power-terms:
544 T b1 = (x * lc) / la;
550 if((lb1 >= tools::log_max_value<T>())
551 || (lb1 <= tools::log_min_value<T>())
552 || (lb2 >= tools::log_max_value<T>())
553 || (lb2 <= tools::log_min_value<T>())
554 || (e1 >= tools::log_max_value<T>())
555 || (e1 <= tools::log_min_value<T>()) )
557 T p = lb1 + lb2 - e1;
564 result *= exp(b * boost::math::log1p(a / lb, pol));
566 result *= pow(b2, b);
569 // and combine the results:
570 result /= sa * sb / sc;
574 *p_derivative = result * pow(y, b);
575 BOOST_ASSERT(*p_derivative >= 0);
580 // Non-normalised, just compute the power:
583 if(result < tools::min_value<T>())
584 return s0; // Safeguard: series can't cope with denorms.
585 ibeta_series_t<T> s(a, b, x, result);
586 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
587 result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
588 policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol);
593 // Continued fraction for the incomplete beta:
596 struct ibeta_fraction2_t
598 typedef std::pair<T, T> result_type;
600 ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {}
602 result_type operator()()
604 T aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x;
605 T denom = (a + 2 * m - 1);
609 bN += (m * (b - m) * x) / (a + 2*m - 1);
610 bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1);
614 return std::make_pair(aN, bN);
622 // Evaluate the incomplete beta via the continued fraction representation:
624 template <class T, class Policy>
625 inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative)
627 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
629 T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
632 *p_derivative = result;
633 BOOST_ASSERT(*p_derivative >= 0);
638 ibeta_fraction2_t<T> f(a, b, x, y);
639 T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>());
640 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
641 BOOST_MATH_INSTRUMENT_VARIABLE(result);
642 return result / fract;
645 // Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x):
647 template <class T, class Policy>
648 T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative)
650 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
652 BOOST_MATH_INSTRUMENT_VARIABLE(k);
654 T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
657 *p_derivative = prefix;
658 BOOST_ASSERT(*p_derivative >= 0);
665 // series summation from 0 to k-1:
666 for(int i = 0; i < k-1; ++i)
668 term *= (a+b+i) * x / (a+i+1);
676 // This function is only needed for the non-regular incomplete beta,
677 // it computes the delta in:
678 // beta(a,b,x) = prefix + delta * beta(a+k,b,x)
679 // it is currently only called for small k.
682 inline T rising_factorial_ratio(T a, T b, int k)
685 // (a)(a+1)(a+2)...(a+k-1)
686 // _______________________
687 // (b)(b+1)(b+2)...(b+k-1)
689 // This is only called with small k, for large k
690 // it is grossly inefficient, do not use outside it's
691 // intended purpose!!!
692 BOOST_MATH_INSTRUMENT_VARIABLE(k);
696 for(int i = 0; i < k; ++i)
697 result *= (a+i) / (b+i);
701 // Routine for a > 15, b < 1
703 // Begin by figuring out how large our table of Pn's should be,
704 // quoted accuracies are "guestimates" based on empiracal observation.
705 // Note that the table size should never exceed the size of our
706 // tables of factorials.
711 // This is likely to be enough for ~35-50 digit accuracy
712 // but it's hard to quantify exactly:
713 BOOST_STATIC_CONSTANT(unsigned, value = 50);
714 BOOST_STATIC_ASSERT(::boost::math::max_factorial<T>::value >= 100);
717 struct Pn_size<float>
719 BOOST_STATIC_CONSTANT(unsigned, value = 15); // ~8-15 digit accuracy
720 BOOST_STATIC_ASSERT(::boost::math::max_factorial<float>::value >= 30);
723 struct Pn_size<double>
725 BOOST_STATIC_CONSTANT(unsigned, value = 30); // 16-20 digit accuracy
726 BOOST_STATIC_ASSERT(::boost::math::max_factorial<double>::value >= 60);
729 struct Pn_size<long double>
731 BOOST_STATIC_CONSTANT(unsigned, value = 50); // ~35-50 digit accuracy
732 BOOST_STATIC_ASSERT(::boost::math::max_factorial<long double>::value >= 100);
735 template <class T, class Policy>
736 T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised)
738 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
741 // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6.
743 // Some values we'll need later, these are Eq 9.1:
749 lx = boost::math::log1p(-y, pol);
753 // and from from 9.2:
755 T h = regularised_gamma_prefix(b, u, pol, lanczos_type());
756 if(h <= tools::min_value<T>())
760 prefix = h / boost::math::tgamma_delta_ratio(a, b, pol);
765 prefix = full_igamma_prefix(b, u, pol) / pow(t, b);
769 // now we need the quantity Pn, unfortunatately this is computed
770 // recursively, and requires a full history of all the previous values
771 // so no choice but to declare a big table and hope it's big enough...
773 T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 }; // see 9.3.
775 // Now an initial value for J, see 9.6:
777 T j = boost::math::gamma_q(b, u, pol) / h;
779 // Now we can start to pull things together and evaluate the sum in Eq 9:
781 T sum = s0 + prefix * j; // Value at N = 0
782 // some variables we'll need:
783 unsigned tnp1 = 1; // 2*N+1
790 for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n)
793 // debugging code, enable this if you want to determine whether
794 // the table of Pn's is large enough...
796 static int max_count = 2;
800 std::cerr << "Max iterations in BGRAT was " << n << std::endl;
804 // begin by evaluating the next Pn from Eq 9.4:
810 for(unsigned m = 1; m < n; ++m)
813 p[n] += mbn * p[n-m] / boost::math::unchecked_factorial<T>(tmp1);
817 p[n] += bm1 / boost::math::unchecked_factorial<T>(tnp1);
819 // Now we want Jn from Jn-1 using Eq 9.6:
821 j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4;
825 // pull it together with Eq 9:
827 T r = prefix * p[n] * j;
831 if(fabs(r) < fabs(tools::epsilon<T>() * sum))
836 if(fabs(r / tools::epsilon<T>()) < fabs(sum))
841 } // template <class T, class Lanczos>T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Lanczos& l, bool normalised)
844 // For integer arguments we can relate the incomplete beta to the
845 // complement of the binomial distribution cdf and use this finite sum.
848 inline T binomial_ccdf(T n, T k, T x, T y)
850 BOOST_MATH_STD_USING // ADL of std names
851 T result = pow(x, n);
853 for(unsigned i = itrunc(T(n - 1)); i > k; --i)
855 term *= ((i + 1) * y) / ((n - i) * x) ;
864 // The incomplete beta function implementation:
865 // This is just a big bunch of spagetti code to divide up the
866 // input range and select the right implementation method for
869 template <class T, class Policy>
870 T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative)
872 static const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)";
873 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
874 BOOST_MATH_STD_USING // for ADL of std math functions.
876 BOOST_MATH_INSTRUMENT_VARIABLE(a);
877 BOOST_MATH_INSTRUMENT_VARIABLE(b);
878 BOOST_MATH_INSTRUMENT_VARIABLE(x);
879 BOOST_MATH_INSTRUMENT_VARIABLE(inv);
880 BOOST_MATH_INSTRUMENT_VARIABLE(normalised);
886 BOOST_ASSERT((p_derivative == 0) || normalised);
889 *p_derivative = -1; // value not set.
891 if((x < 0) || (x > 1))
892 policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
897 policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
899 policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
900 // extend to a few very special cases:
904 policies::raise_domain_error<T>(function, "The arguments a and b to the incomplete beta function cannot both be zero, with x=%1%.", x, pol);
917 policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
919 policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
926 *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
928 return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0));
934 *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
936 return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
938 if((a == 0.5f) && (b == 0.5f))
940 // We have an arcsine distribution:
943 *p_derivative = (invert ? -1 : 1) / constants::pi<T>() * sqrt(y * x);
945 T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>();
947 p *= constants::pi<T>();
959 // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
964 *p_derivative = invert ? -1 : 1;
965 return invert ? y : x;
970 *p_derivative = (invert ? -1 : 1) * a * pow(x, a - 1);
974 p = invert ? T(-expm1(a * log1p(-y))) : T(exp(a * log1p(-y)));
976 p = invert ? T(-powm1(x, a)) : T(pow(x, a));
982 if((std::min)(a, b) <= 1)
989 BOOST_MATH_INSTRUMENT_VARIABLE(invert);
991 if((std::max)(a, b) <= 1)
994 if((a >= (std::min)(T(0.2), b)) || (pow(x, a) <= 0.9))
998 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
999 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1003 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1005 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1006 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1018 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1019 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1023 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1025 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1026 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1031 // Sidestep on a, and then use the series representation:
1035 prefix = rising_factorial_ratio(T(a+b), a, 20);
1041 fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
1044 fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1045 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1049 fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1051 fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1052 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1059 // One of a, b < 1 only:
1060 if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7)))
1064 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1065 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1069 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1071 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1072 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1085 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1086 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1090 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1092 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1093 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1100 fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised);
1101 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1105 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1107 fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised);
1108 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1113 // Sidestep to improve errors:
1117 prefix = rising_factorial_ratio(T(a+b), a, 20);
1123 fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
1124 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1127 fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1128 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1132 fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1134 fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1135 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1147 lambda = a - (a + b) * x;
1151 lambda = (a + b) * y - b;
1158 BOOST_MATH_INSTRUMENT_VARIABLE(invert);
1163 if((floor(a) == a) && (floor(b) == b) && (a < (std::numeric_limits<int>::max)() - 100))
1165 // relate to the binomial distribution and use a finite sum:
1168 fract = binomial_ccdf(n, k, x, y);
1170 fract *= boost::math::beta(a, b, pol);
1171 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1173 else if(b * x <= 0.7)
1177 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1178 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1182 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1184 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1185 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1190 // sidestep so we can use the series representation:
1191 int n = itrunc(T(floor(b)), pol);
1198 prefix = rising_factorial_ratio(T(a+bbar), bbar, n);
1204 fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
1205 fract = beta_small_b_large_a_series(a, bbar, x, y, fract, T(1), pol, normalised);
1207 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1211 // the formula here for the non-normalised case is tricky to figure
1212 // out (for me!!), and requires two pochhammer calculations rather
1213 // than one, so leave it for now....
1214 int n = itrunc(T(floor(b)), pol);
1221 fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
1222 fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(0));
1224 fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1225 //fract = ibeta_series(a+20, bbar, x, fract, l, normalised, p_derivative, y);
1226 fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised);
1232 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1236 fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
1237 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1242 fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
1243 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1248 if(*p_derivative < 0)
1250 *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol);
1254 if(*p_derivative != 0)
1256 if((tools::max_value<T>() * div < *p_derivative))
1258 // overflow, return an arbitarily large value:
1259 *p_derivative = tools::max_value<T>() / 2;
1263 *p_derivative /= div;
1267 return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
1268 } // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised)
1270 template <class T, class Policy>
1271 inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised)
1273 return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(0));
1276 template <class T, class Policy>
1277 T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
1279 static const char* function = "ibeta_derivative<%1%>(%1%,%1%,%1%)";
1281 // start with the usual error checks:
1284 policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
1286 policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
1287 if((x < 0) || (x > 1))
1288 policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
1290 // Now the corner cases:
1294 return (a > 1) ? 0 :
1295 (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol);
1299 return (b > 1) ? 0 :
1300 (b == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol);
1303 // Now the regular cases:
1305 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1306 T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol);
1312 if((tools::max_value<T>() * y < f1))
1315 return policies::raise_overflow_error<T>(function, 0, pol);
1323 // Some forwarding functions that dis-ambiguate the third argument type:
1325 template <class RT1, class RT2, class Policy>
1326 inline typename tools::promote_args<RT1, RT2>::type
1327 beta(RT1 a, RT2 b, const Policy&, const mpl::true_*)
1329 BOOST_FPU_EXCEPTION_GUARD
1330 typedef typename tools::promote_args<RT1, RT2>::type result_type;
1331 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1332 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1333 typedef typename policies::normalise<
1335 policies::promote_float<false>,
1336 policies::promote_double<false>,
1337 policies::discrete_quantile<>,
1338 policies::assert_undefined<> >::type forwarding_policy;
1340 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::beta_imp(static_cast<value_type>(a), static_cast<value_type>(b), evaluation_type(), forwarding_policy()), "boost::math::beta<%1%>(%1%,%1%)");
1342 template <class RT1, class RT2, class RT3>
1343 inline typename tools::promote_args<RT1, RT2, RT3>::type
1344 beta(RT1 a, RT2 b, RT3 x, const mpl::false_*)
1346 return boost::math::beta(a, b, x, policies::policy<>());
1348 } // namespace detail
1351 // The actual function entry-points now follow, these just figure out
1352 // which Lanczos approximation to use
1353 // and forward to the implementation functions:
1355 template <class RT1, class RT2, class A>
1356 inline typename tools::promote_args<RT1, RT2, A>::type
1357 beta(RT1 a, RT2 b, A arg)
1359 typedef typename policies::is_policy<A>::type tag;
1360 return boost::math::detail::beta(a, b, arg, static_cast<tag*>(0));
1363 template <class RT1, class RT2>
1364 inline typename tools::promote_args<RT1, RT2>::type
1367 return boost::math::beta(a, b, policies::policy<>());
1370 template <class RT1, class RT2, class RT3, class Policy>
1371 inline typename tools::promote_args<RT1, RT2, RT3>::type
1372 beta(RT1 a, RT2 b, RT3 x, const Policy&)
1374 BOOST_FPU_EXCEPTION_GUARD
1375 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1376 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1377 typedef typename policies::normalise<
1379 policies::promote_float<false>,
1380 policies::promote_double<false>,
1381 policies::discrete_quantile<>,
1382 policies::assert_undefined<> >::type forwarding_policy;
1384 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, false), "boost::math::beta<%1%>(%1%,%1%,%1%)");
1387 template <class RT1, class RT2, class RT3, class Policy>
1388 inline typename tools::promote_args<RT1, RT2, RT3>::type
1389 betac(RT1 a, RT2 b, RT3 x, const Policy&)
1391 BOOST_FPU_EXCEPTION_GUARD
1392 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1393 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1394 typedef typename policies::normalise<
1396 policies::promote_float<false>,
1397 policies::promote_double<false>,
1398 policies::discrete_quantile<>,
1399 policies::assert_undefined<> >::type forwarding_policy;
1401 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, false), "boost::math::betac<%1%>(%1%,%1%,%1%)");
1403 template <class RT1, class RT2, class RT3>
1404 inline typename tools::promote_args<RT1, RT2, RT3>::type
1405 betac(RT1 a, RT2 b, RT3 x)
1407 return boost::math::betac(a, b, x, policies::policy<>());
1410 template <class RT1, class RT2, class RT3, class Policy>
1411 inline typename tools::promote_args<RT1, RT2, RT3>::type
1412 ibeta(RT1 a, RT2 b, RT3 x, const Policy&)
1414 BOOST_FPU_EXCEPTION_GUARD
1415 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1416 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1417 typedef typename policies::normalise<
1419 policies::promote_float<false>,
1420 policies::promote_double<false>,
1421 policies::discrete_quantile<>,
1422 policies::assert_undefined<> >::type forwarding_policy;
1424 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, true), "boost::math::ibeta<%1%>(%1%,%1%,%1%)");
1426 template <class RT1, class RT2, class RT3>
1427 inline typename tools::promote_args<RT1, RT2, RT3>::type
1428 ibeta(RT1 a, RT2 b, RT3 x)
1430 return boost::math::ibeta(a, b, x, policies::policy<>());
1433 template <class RT1, class RT2, class RT3, class Policy>
1434 inline typename tools::promote_args<RT1, RT2, RT3>::type
1435 ibetac(RT1 a, RT2 b, RT3 x, const Policy&)
1437 BOOST_FPU_EXCEPTION_GUARD
1438 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1439 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1440 typedef typename policies::normalise<
1442 policies::promote_float<false>,
1443 policies::promote_double<false>,
1444 policies::discrete_quantile<>,
1445 policies::assert_undefined<> >::type forwarding_policy;
1447 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, true), "boost::math::ibetac<%1%>(%1%,%1%,%1%)");
1449 template <class RT1, class RT2, class RT3>
1450 inline typename tools::promote_args<RT1, RT2, RT3>::type
1451 ibetac(RT1 a, RT2 b, RT3 x)
1453 return boost::math::ibetac(a, b, x, policies::policy<>());
1456 template <class RT1, class RT2, class RT3, class Policy>
1457 inline typename tools::promote_args<RT1, RT2, RT3>::type
1458 ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy&)
1460 BOOST_FPU_EXCEPTION_GUARD
1461 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1462 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1463 typedef typename policies::normalise<
1465 policies::promote_float<false>,
1466 policies::promote_double<false>,
1467 policies::discrete_quantile<>,
1468 policies::assert_undefined<> >::type forwarding_policy;
1470 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy()), "boost::math::ibeta_derivative<%1%>(%1%,%1%,%1%)");
1472 template <class RT1, class RT2, class RT3>
1473 inline typename tools::promote_args<RT1, RT2, RT3>::type
1474 ibeta_derivative(RT1 a, RT2 b, RT3 x)
1476 return boost::math::ibeta_derivative(a, b, x, policies::policy<>());
1480 } // namespace boost
1482 #include <boost/math/special_functions/detail/ibeta_inverse.hpp>
1483 #include <boost/math/special_functions/detail/ibeta_inv_ab.hpp>
1485 #endif // BOOST_MATH_SPECIAL_BETA_HPP