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/binomial.hpp>
17 #include <boost/math/special_functions/factorials.hpp>
18 #include <boost/math/special_functions/erf.hpp>
19 #include <boost/math/special_functions/log1p.hpp>
20 #include <boost/math/special_functions/expm1.hpp>
21 #include <boost/math/special_functions/trunc.hpp>
22 #include <boost/math/tools/roots.hpp>
23 #include <boost/static_assert.hpp>
24 #include <boost/config/no_tr1/cmath.hpp>
26 namespace boost{ namespace math{
31 // Implementation of Beta(a,b) using the Lanczos approximation:
33 template <class T, class Lanczos, class Policy>
34 T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
36 BOOST_MATH_STD_USING // for ADL of std names
39 return 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);
41 return 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);
49 if((c == a) && (b < tools::epsilon<T>()))
51 else if((c == b) && (a < tools::epsilon<T>()))
57 else if(c < tools::epsilon<T>())
66 // This code appears to be no longer necessary: it was
67 // used to offset errors introduced from the Lanczos
68 // approximation, but the current Lanczos approximations
69 // are sufficiently accurate for all z that we can ditch
70 // this. It remains in the file for future reference...
72 // If a or b are less than 1, shift to greater than 1:
90 // Lanczos calculation:
91 T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
92 T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
93 T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
94 result = Lanczos::lanczos_sum_expG_scaled(a) * (Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c));
95 T ambh = a - 0.5f - b;
96 if((fabs(b * ambh) < (cgh * 100)) && (a > 100))
98 // Special case where the base of the power term is close to 1
99 // compute (1+x)^y instead:
100 result *= exp(ambh * boost::math::log1p(-b / cgh, pol));
104 result *= pow(agh / cgh, a - T(0.5) - b);
107 // this avoids possible overflow, but appears to be marginally less accurate:
108 result *= pow((agh / cgh) * (bgh / cgh), b);
110 result *= pow((agh * bgh) / (cgh * cgh), b);
111 result *= sqrt(boost::math::constants::e<T>() / bgh);
113 // If a and b were originally less than 1 we need to scale the result:
117 } // template <class T, class Lanczos> beta_imp(T a, T b, const Lanczos&)
120 // Generic implementation of Beta(a,b) without Lanczos approximation support
121 // (Caution this is slow!!!):
123 template <class T, class Policy>
124 T beta_imp(T a, T b, const lanczos::undefined_lanczos& /* l */, const Policy& pol)
129 return 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);
131 return 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);
139 if((c == a) && (b < tools::epsilon<T>()))
140 return boost::math::tgamma(b, pol);
141 else if((c == b) && (a < tools::epsilon<T>()))
142 return boost::math::tgamma(a, pol);
148 // shift to a and b > 1 if required:
164 // set integration limits:
165 T la = (std::max)(T(10), a);
166 T lb = (std::max)(T(10), b);
167 T lc = (std::max)(T(10), T(a+b));
169 // calculate the fraction parts:
170 T sa = detail::lower_gamma_series(a, la, pol) / a;
171 sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
172 T sb = detail::lower_gamma_series(b, lb, pol) / b;
173 sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
174 T sc = detail::lower_gamma_series(c, lc, pol) / c;
175 sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
177 // and the exponent part:
178 result = exp(lc - la - lb) * pow(la/lc, a) * pow(lb/lc, b);
181 result *= sa * sb / sc;
183 // if a and b were originally less than 1 we need to scale the result:
187 } // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l)
191 // Compute the leading power terms in the incomplete Beta:
193 // (x^a)(y^b)/Beta(a,b) when normalised, and
194 // (x^a)(y^b) otherwise.
196 // Almost all of the error in the incomplete beta comes from this
197 // function: particularly when a and b are large. Computing large
198 // powers are *hard* though, and using logarithms just leads to
199 // horrendous cancellation errors.
201 template <class T, class Lanczos, class Policy>
202 T ibeta_power_terms(T a,
210 const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
216 // can we do better here?
217 return pow(x, a) * pow(y, b);
224 // combine power terms with Lanczos approximation:
225 T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
226 T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
227 T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
228 result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
230 // combine with the leftover terms from the Lanczos approximation:
231 result *= sqrt(bgh / boost::math::constants::e<T>());
232 result *= sqrt(agh / cgh);
234 // l1 and l2 are the base of the exponents minus one:
235 T l1 = (x * b - y * agh) / agh;
236 T l2 = (y * a - x * bgh) / bgh;
237 if(((std::min)(fabs(l1), fabs(l2)) < 0.2))
239 // when the base of the exponent is very near 1 we get really
240 // gross errors unless extra care is taken:
241 if((l1 * l2 > 0) || ((std::min)(a, b) < 1))
244 // This first branch handles the simple cases where either:
246 // * The two power terms both go in the same direction
247 // (towards zero or towards infinity). In this case if either
248 // term overflows or underflows, then the product of the two must
250 // *Alternatively if one exponent is less than one, then we
251 // can't productively use it to eliminate overflow or underflow
252 // from the other term. Problems with spurious overflow/underflow
253 // can't be ruled out in this case, but it is *very* unlikely
254 // since one of the power terms will evaluate to a number close to 1.
258 result *= exp(a * boost::math::log1p(l1, pol));
259 BOOST_MATH_INSTRUMENT_VARIABLE(result);
263 result *= pow((x * cgh) / agh, a);
264 BOOST_MATH_INSTRUMENT_VARIABLE(result);
268 result *= exp(b * boost::math::log1p(l2, pol));
269 BOOST_MATH_INSTRUMENT_VARIABLE(result);
273 result *= pow((y * cgh) / bgh, b);
274 BOOST_MATH_INSTRUMENT_VARIABLE(result);
277 else if((std::max)(fabs(l1), fabs(l2)) < 0.5)
280 // Both exponents are near one and both the exponents are
281 // greater than one and further these two
282 // power terms tend in opposite directions (one towards zero,
283 // the other towards infinity), so we have to combine the terms
284 // to avoid any risk of overflow or underflow.
286 // We do this by moving one power term inside the other, we have:
288 // (1 + l1)^a * (1 + l2)^b
289 // = ((1 + l1)*(1 + l2)^(b/a))^a
290 // = (1 + l1 + l3 + l1*l3)^a ; l3 = (1 + l2)^(b/a) - 1
291 // = exp((b/a) * log(1 + l2)) - 1
293 // The tricky bit is deciding which term to move inside :-)
294 // By preference we move the larger term inside, so that the
295 // size of the largest exponent is reduced. However, that can
296 // only be done as long as l3 (see above) is also small.
298 bool small_a = a < b;
300 if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1)))
302 T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol);
303 l3 = l1 + l3 + l3 * l1;
304 l3 = a * boost::math::log1p(l3, pol);
306 BOOST_MATH_INSTRUMENT_VARIABLE(result);
310 T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol);
311 l3 = l2 + l3 + l3 * l2;
312 l3 = b * boost::math::log1p(l3, pol);
314 BOOST_MATH_INSTRUMENT_VARIABLE(result);
317 else if(fabs(l1) < fabs(l2))
319 // First base near 1 only:
320 T l = a * boost::math::log1p(l1, pol)
321 + b * log((y * cgh) / bgh);
322 if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>()))
325 if(l >= tools::log_max_value<T>())
326 return policies::raise_overflow_error<T>(function, 0, pol);
331 BOOST_MATH_INSTRUMENT_VARIABLE(result);
335 // Second base near 1 only:
336 T l = b * boost::math::log1p(l2, pol)
337 + a * log((x * cgh) / agh);
338 if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>()))
341 if(l >= tools::log_max_value<T>())
342 return policies::raise_overflow_error<T>(function, 0, pol);
347 BOOST_MATH_INSTRUMENT_VARIABLE(result);
353 T b1 = (x * cgh) / agh;
354 T b2 = (y * cgh) / bgh;
357 BOOST_MATH_INSTRUMENT_VARIABLE(b1);
358 BOOST_MATH_INSTRUMENT_VARIABLE(b2);
359 BOOST_MATH_INSTRUMENT_VARIABLE(l1);
360 BOOST_MATH_INSTRUMENT_VARIABLE(l2);
361 if((l1 >= tools::log_max_value<T>())
362 || (l1 <= tools::log_min_value<T>())
363 || (l2 >= tools::log_max_value<T>())
364 || (l2 <= tools::log_min_value<T>())
367 // Oops, under/overflow, sidestep if we can:
370 T p1 = pow(b2, b / a);
371 T l3 = a * (log(b1) + log(p1));
372 if((l3 < tools::log_max_value<T>())
373 && (l3 > tools::log_min_value<T>()))
375 result *= pow(p1 * b1, a);
379 l2 += l1 + log(result);
380 if(l2 >= tools::log_max_value<T>())
381 return policies::raise_overflow_error<T>(function, 0, pol);
387 T p1 = pow(b1, a / b);
388 T l3 = (log(p1) + log(b2)) * b;
389 if((l3 < tools::log_max_value<T>())
390 && (l3 > tools::log_min_value<T>()))
392 result *= pow(p1 * b2, b);
396 l2 += l1 + log(result);
397 if(l2 >= tools::log_max_value<T>())
398 return policies::raise_overflow_error<T>(function, 0, pol);
402 BOOST_MATH_INSTRUMENT_VARIABLE(result);
406 // finally the normal case:
407 result *= pow(b1, a) * pow(b2, b);
408 BOOST_MATH_INSTRUMENT_VARIABLE(result);
412 BOOST_MATH_INSTRUMENT_VARIABLE(result);
417 // Compute the leading power terms in the incomplete Beta:
419 // (x^a)(y^b)/Beta(a,b) when normalised, and
420 // (x^a)(y^b) otherwise.
422 // Almost all of the error in the incomplete beta comes from this
423 // function: particularly when a and b are large. Computing large
424 // powers are *hard* though, and using logarithms just leads to
425 // horrendous cancellation errors.
427 // This version is generic, slow, and does not use the Lanczos approximation.
429 template <class T, class Policy>
430 T ibeta_power_terms(T a,
434 const boost::math::lanczos::undefined_lanczos&,
438 const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
444 return pow(x, a) * pow(y, b);
447 T result= 0; // assignment here silences warnings later
451 // integration limits for the gamma functions:
452 //T la = (std::max)(T(10), a);
453 //T lb = (std::max)(T(10), b);
454 //T lc = (std::max)(T(10), a+b);
458 // gamma function partials:
459 T sa = detail::lower_gamma_series(a, la, pol) / a;
460 sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
461 T sb = detail::lower_gamma_series(b, lb, pol) / b;
462 sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
463 T sc = detail::lower_gamma_series(c, lc, pol) / c;
464 sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
465 // gamma function powers combined with incomplete beta powers:
467 T b1 = (x * lc) / la;
468 T b2 = (y * lc) / lb;
469 T e1 = -5; // lc - la - lb;
473 if((lb1 >= tools::log_max_value<T>())
474 || (lb1 <= tools::log_min_value<T>())
475 || (lb2 >= tools::log_max_value<T>())
476 || (lb2 <= tools::log_min_value<T>())
477 || (e1 >= tools::log_max_value<T>())
478 || (e1 <= tools::log_min_value<T>())
481 result = exp(lb1 + lb2 - e1 + log(prefix));
486 p1 = (x * b - y * la) / la;
488 p1 = exp(a * boost::math::log1p(p1, pol));
491 p2 = (y * a - x * lb) / lb;
493 p2 = exp(b * boost::math::log1p(p2, pol));
497 result = prefix * p1 * (p2 / p3);
499 // and combine with the remaining gamma function components:
500 result /= sa * sb / sc;
505 // Series approximation to the incomplete beta:
508 struct ibeta_series_t
510 typedef T result_type;
511 ibeta_series_t(T a_, T b_, T x_, T mult) : result(mult), x(x_), apn(a_), poch(1-b_), n(1) {}
516 result *= poch * x / n;
522 T result, x, apn, poch;
526 template <class T, class Lanczos, class Policy>
527 T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol)
533 BOOST_ASSERT((p_derivative == 0) || normalised);
539 // incomplete beta power term, combined with the Lanczos approximation:
540 T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
541 T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
542 T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
543 result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
545 T l1 = log(cgh / bgh) * (b - 0.5f);
546 T l2 = log(x * cgh / agh) * a;
548 // Check for over/underflow in the power terms:
550 if((l1 > tools::log_min_value<T>())
551 && (l1 < tools::log_max_value<T>())
552 && (l2 > tools::log_min_value<T>())
553 && (l2 < tools::log_max_value<T>()))
556 result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol));
558 result *= pow(cgh / bgh, b - 0.5f);
559 result *= pow(x * cgh / agh, a);
560 result *= sqrt(agh / boost::math::constants::e<T>());
564 *p_derivative = result * pow(y, b);
565 BOOST_ASSERT(*p_derivative >= 0);
571 // Oh dear, we need logs, and this *will* cancel:
573 result = log(result) + l1 + l2 + (log(agh) - 1) / 2;
575 *p_derivative = exp(result + b * log(y));
576 result = exp(result);
581 // Non-normalised, just compute the power:
584 if(result < tools::min_value<T>())
585 return s0; // Safeguard: series can't cope with denorms.
586 ibeta_series_t<T> s(a, b, x, result);
587 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
588 result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
589 policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol);
593 // Incomplete Beta series again, this time without Lanczos support:
595 template <class T, class Policy>
596 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)
601 BOOST_ASSERT((p_derivative == 0) || normalised);
607 // figure out integration limits for the gamma function:
608 //T la = (std::max)(T(10), a);
609 //T lb = (std::max)(T(10), b);
610 //T lc = (std::max)(T(10), a+b);
615 // calculate the gamma parts:
616 T sa = detail::lower_gamma_series(a, la, pol) / a;
617 sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
618 T sb = detail::lower_gamma_series(b, lb, pol) / b;
619 sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
620 T sc = detail::lower_gamma_series(c, lc, pol) / c;
621 sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
623 // and their combined power-terms:
624 T b1 = (x * lc) / la;
630 if((lb1 >= tools::log_max_value<T>())
631 || (lb1 <= tools::log_min_value<T>())
632 || (lb2 >= tools::log_max_value<T>())
633 || (lb2 <= tools::log_min_value<T>())
634 || (e1 >= tools::log_max_value<T>())
635 || (e1 <= tools::log_min_value<T>()) )
637 T p = lb1 + lb2 - e1;
644 result *= exp(b * boost::math::log1p(a / lb, pol));
646 result *= pow(b2, b);
649 // and combine the results:
650 result /= sa * sb / sc;
654 *p_derivative = result * pow(y, b);
655 BOOST_ASSERT(*p_derivative >= 0);
660 // Non-normalised, just compute the power:
663 if(result < tools::min_value<T>())
664 return s0; // Safeguard: series can't cope with denorms.
665 ibeta_series_t<T> s(a, b, x, result);
666 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
667 result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
668 policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol);
673 // Continued fraction for the incomplete beta:
676 struct ibeta_fraction2_t
678 typedef std::pair<T, T> result_type;
680 ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {}
682 result_type operator()()
684 T aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x;
685 T denom = (a + 2 * m - 1);
688 T bN = static_cast<T>(m);
689 bN += (m * (b - m) * x) / (a + 2*m - 1);
690 bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1);
694 return std::make_pair(aN, bN);
702 // Evaluate the incomplete beta via the continued fraction representation:
704 template <class T, class Policy>
705 inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative)
707 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
709 T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
712 *p_derivative = result;
713 BOOST_ASSERT(*p_derivative >= 0);
718 ibeta_fraction2_t<T> f(a, b, x, y);
719 T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>());
720 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
721 BOOST_MATH_INSTRUMENT_VARIABLE(result);
722 return result / fract;
725 // Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x):
727 template <class T, class Policy>
728 T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative)
730 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
732 BOOST_MATH_INSTRUMENT_VARIABLE(k);
734 T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
737 *p_derivative = prefix;
738 BOOST_ASSERT(*p_derivative >= 0);
745 // series summation from 0 to k-1:
746 for(int i = 0; i < k-1; ++i)
748 term *= (a+b+i) * x / (a+i+1);
756 // This function is only needed for the non-regular incomplete beta,
757 // it computes the delta in:
758 // beta(a,b,x) = prefix + delta * beta(a+k,b,x)
759 // it is currently only called for small k.
762 inline T rising_factorial_ratio(T a, T b, int k)
765 // (a)(a+1)(a+2)...(a+k-1)
766 // _______________________
767 // (b)(b+1)(b+2)...(b+k-1)
769 // This is only called with small k, for large k
770 // it is grossly inefficient, do not use outside it's
771 // intended purpose!!!
772 BOOST_MATH_INSTRUMENT_VARIABLE(k);
776 for(int i = 0; i < k; ++i)
777 result *= (a+i) / (b+i);
781 // Routine for a > 15, b < 1
783 // Begin by figuring out how large our table of Pn's should be,
784 // quoted accuracies are "guestimates" based on empiracal observation.
785 // Note that the table size should never exceed the size of our
786 // tables of factorials.
791 // This is likely to be enough for ~35-50 digit accuracy
792 // but it's hard to quantify exactly:
793 BOOST_STATIC_CONSTANT(unsigned, value = 50);
794 BOOST_STATIC_ASSERT(::boost::math::max_factorial<T>::value >= 100);
797 struct Pn_size<float>
799 BOOST_STATIC_CONSTANT(unsigned, value = 15); // ~8-15 digit accuracy
800 BOOST_STATIC_ASSERT(::boost::math::max_factorial<float>::value >= 30);
803 struct Pn_size<double>
805 BOOST_STATIC_CONSTANT(unsigned, value = 30); // 16-20 digit accuracy
806 BOOST_STATIC_ASSERT(::boost::math::max_factorial<double>::value >= 60);
809 struct Pn_size<long double>
811 BOOST_STATIC_CONSTANT(unsigned, value = 50); // ~35-50 digit accuracy
812 BOOST_STATIC_ASSERT(::boost::math::max_factorial<long double>::value >= 100);
815 template <class T, class Policy>
816 T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised)
818 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
821 // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6.
823 // Some values we'll need later, these are Eq 9.1:
829 lx = boost::math::log1p(-y, pol);
833 // and from from 9.2:
835 T h = regularised_gamma_prefix(b, u, pol, lanczos_type());
836 if(h <= tools::min_value<T>())
840 prefix = h / boost::math::tgamma_delta_ratio(a, b, pol);
845 prefix = full_igamma_prefix(b, u, pol) / pow(t, b);
849 // now we need the quantity Pn, unfortunatately this is computed
850 // recursively, and requires a full history of all the previous values
851 // so no choice but to declare a big table and hope it's big enough...
853 T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 }; // see 9.3.
855 // Now an initial value for J, see 9.6:
857 T j = boost::math::gamma_q(b, u, pol) / h;
859 // Now we can start to pull things together and evaluate the sum in Eq 9:
861 T sum = s0 + prefix * j; // Value at N = 0
862 // some variables we'll need:
863 unsigned tnp1 = 1; // 2*N+1
870 for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n)
873 // debugging code, enable this if you want to determine whether
874 // the table of Pn's is large enough...
876 static int max_count = 2;
880 std::cerr << "Max iterations in BGRAT was " << n << std::endl;
884 // begin by evaluating the next Pn from Eq 9.4:
890 for(unsigned m = 1; m < n; ++m)
893 p[n] += mbn * p[n-m] / boost::math::unchecked_factorial<T>(tmp1);
897 p[n] += bm1 / boost::math::unchecked_factorial<T>(tnp1);
899 // Now we want Jn from Jn-1 using Eq 9.6:
901 j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4;
905 // pull it together with Eq 9:
907 T r = prefix * p[n] * j;
911 if(fabs(r) < fabs(tools::epsilon<T>() * sum))
916 if(fabs(r / tools::epsilon<T>()) < fabs(sum))
921 } // 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)
924 // For integer arguments we can relate the incomplete beta to the
925 // complement of the binomial distribution cdf and use this finite sum.
928 T binomial_ccdf(T n, T k, T x, T y)
930 BOOST_MATH_STD_USING // ADL of std names
932 T result = pow(x, n);
934 if(result > tools::min_value<T>())
937 for(unsigned i = itrunc(T(n - 1)); i > k; --i)
939 term *= ((i + 1) * y) / ((n - i) * x);
945 // First term underflows so we need to start at the mode of the
946 // distribution and work outwards:
947 int start = itrunc(n * x);
949 start = itrunc(k + 2);
950 result = pow(x, start) * pow(y, n - start) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start));
953 // OK, starting slightly above the mode didn't work,
954 // we'll have to sum the terms the old fashioned way:
955 for(unsigned i = start - 1; i > k; --i)
957 result += pow(x, (int)i) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i));
963 T start_term = result;
964 for(unsigned i = start - 1; i > k; --i)
966 term *= ((i + 1) * y) / ((n - i) * x);
970 for(unsigned i = start + 1; i <= n; ++i)
972 term *= (n - i + 1) * x / (i * y);
983 // The incomplete beta function implementation:
984 // This is just a big bunch of spagetti code to divide up the
985 // input range and select the right implementation method for
988 template <class T, class Policy>
989 T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative)
991 static const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)";
992 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
993 BOOST_MATH_STD_USING // for ADL of std math functions.
995 BOOST_MATH_INSTRUMENT_VARIABLE(a);
996 BOOST_MATH_INSTRUMENT_VARIABLE(b);
997 BOOST_MATH_INSTRUMENT_VARIABLE(x);
998 BOOST_MATH_INSTRUMENT_VARIABLE(inv);
999 BOOST_MATH_INSTRUMENT_VARIABLE(normalised);
1005 BOOST_ASSERT((p_derivative == 0) || normalised);
1008 *p_derivative = -1; // value not set.
1010 if((x < 0) || (x > 1))
1011 return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
1016 return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
1018 return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
1019 // extend to a few very special cases:
1023 return 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);
1025 return static_cast<T>(inv ? 0 : 1);
1030 return static_cast<T>(inv ? 1 : 0);
1036 return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
1038 return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
1045 *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
1047 return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0));
1053 *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
1055 return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
1057 if((a == 0.5f) && (b == 0.5f))
1059 // We have an arcsine distribution:
1062 *p_derivative = 1 / constants::pi<T>() * sqrt(y * x);
1064 T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>();
1066 p *= constants::pi<T>();
1078 // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
1084 return invert ? y : x;
1089 *p_derivative = a * pow(x, a - 1);
1093 p = invert ? T(-boost::math::expm1(a * boost::math::log1p(-y, pol), pol)) : T(exp(a * boost::math::log1p(-y, pol)));
1095 p = invert ? T(-boost::math::powm1(x, a, pol)) : T(pow(x, a));
1101 if((std::min)(a, b) <= 1)
1108 BOOST_MATH_INSTRUMENT_VARIABLE(invert);
1110 if((std::max)(a, b) <= 1)
1113 if((a >= (std::min)(T(0.2), b)) || (pow(x, a) <= 0.9))
1117 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1118 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1122 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1124 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1125 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1137 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1138 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1142 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1144 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1145 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1150 // Sidestep on a, and then use the series representation:
1154 prefix = rising_factorial_ratio(T(a+b), a, 20);
1160 fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
1163 fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1164 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1168 fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1170 fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1171 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1178 // One of a, b < 1 only:
1179 if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7)))
1183 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1184 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1188 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1190 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1191 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1204 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1205 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1209 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1211 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1212 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1219 fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised);
1220 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1224 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1226 fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised);
1227 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1232 // Sidestep to improve errors:
1236 prefix = rising_factorial_ratio(T(a+b), a, 20);
1242 fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
1243 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1246 fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1247 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1251 fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1253 fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1254 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1266 lambda = a - (a + b) * x;
1270 lambda = (a + b) * y - b;
1277 BOOST_MATH_INSTRUMENT_VARIABLE(invert);
1282 if((floor(a) == a) && (floor(b) == b) && (a < (std::numeric_limits<int>::max)() - 100) && (y != 1))
1284 // relate to the binomial distribution and use a finite sum:
1287 fract = binomial_ccdf(n, k, x, y);
1289 fract *= boost::math::beta(a, b, pol);
1290 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1292 else if(b * x <= 0.7)
1296 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1297 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1301 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1303 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1304 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1309 // sidestep so we can use the series representation:
1310 int n = itrunc(T(floor(b)), pol);
1317 prefix = rising_factorial_ratio(T(a+bbar), bbar, n);
1323 fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
1324 fract = beta_small_b_large_a_series(a, bbar, x, y, fract, T(1), pol, normalised);
1326 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1330 // The formula here for the non-normalised case is tricky to figure
1331 // out (for me!!), and requires two pochhammer calculations rather
1332 // than one, so leave it for now and only use this in the normalized case....
1333 int n = itrunc(T(floor(b)), pol);
1340 fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
1341 fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(0));
1343 fract -= 1; // Note this line would need changing if we ever enable this branch in non-normalized case
1344 fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised);
1350 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1354 fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
1355 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1360 fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
1361 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1366 if(*p_derivative < 0)
1368 *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol);
1372 if(*p_derivative != 0)
1374 if((tools::max_value<T>() * div < *p_derivative))
1376 // overflow, return an arbitarily large value:
1377 *p_derivative = tools::max_value<T>() / 2;
1381 *p_derivative /= div;
1385 return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
1386 } // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised)
1388 template <class T, class Policy>
1389 inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised)
1391 return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(0));
1394 template <class T, class Policy>
1395 T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
1397 static const char* function = "ibeta_derivative<%1%>(%1%,%1%,%1%)";
1399 // start with the usual error checks:
1402 return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
1404 return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
1405 if((x < 0) || (x > 1))
1406 return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
1408 // Now the corner cases:
1412 return (a > 1) ? 0 :
1413 (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol);
1417 return (b > 1) ? 0 :
1418 (b == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol);
1421 // Now the regular cases:
1423 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1425 T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function);
1429 // Some forwarding functions that dis-ambiguate the third argument type:
1431 template <class RT1, class RT2, class Policy>
1432 inline typename tools::promote_args<RT1, RT2>::type
1433 beta(RT1 a, RT2 b, const Policy&, const mpl::true_*)
1435 BOOST_FPU_EXCEPTION_GUARD
1436 typedef typename tools::promote_args<RT1, RT2>::type result_type;
1437 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1438 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1439 typedef typename policies::normalise<
1441 policies::promote_float<false>,
1442 policies::promote_double<false>,
1443 policies::discrete_quantile<>,
1444 policies::assert_undefined<> >::type forwarding_policy;
1446 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%)");
1448 template <class RT1, class RT2, class RT3>
1449 inline typename tools::promote_args<RT1, RT2, RT3>::type
1450 beta(RT1 a, RT2 b, RT3 x, const mpl::false_*)
1452 return boost::math::beta(a, b, x, policies::policy<>());
1454 } // namespace detail
1457 // The actual function entry-points now follow, these just figure out
1458 // which Lanczos approximation to use
1459 // and forward to the implementation functions:
1461 template <class RT1, class RT2, class A>
1462 inline typename tools::promote_args<RT1, RT2, A>::type
1463 beta(RT1 a, RT2 b, A arg)
1465 typedef typename policies::is_policy<A>::type tag;
1466 return boost::math::detail::beta(a, b, arg, static_cast<tag*>(0));
1469 template <class RT1, class RT2>
1470 inline typename tools::promote_args<RT1, RT2>::type
1473 return boost::math::beta(a, b, policies::policy<>());
1476 template <class RT1, class RT2, class RT3, class Policy>
1477 inline typename tools::promote_args<RT1, RT2, RT3>::type
1478 beta(RT1 a, RT2 b, RT3 x, const Policy&)
1480 BOOST_FPU_EXCEPTION_GUARD
1481 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1482 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1483 typedef typename policies::normalise<
1485 policies::promote_float<false>,
1486 policies::promote_double<false>,
1487 policies::discrete_quantile<>,
1488 policies::assert_undefined<> >::type forwarding_policy;
1490 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%)");
1493 template <class RT1, class RT2, class RT3, class Policy>
1494 inline typename tools::promote_args<RT1, RT2, RT3>::type
1495 betac(RT1 a, RT2 b, RT3 x, const Policy&)
1497 BOOST_FPU_EXCEPTION_GUARD
1498 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1499 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1500 typedef typename policies::normalise<
1502 policies::promote_float<false>,
1503 policies::promote_double<false>,
1504 policies::discrete_quantile<>,
1505 policies::assert_undefined<> >::type forwarding_policy;
1507 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%)");
1509 template <class RT1, class RT2, class RT3>
1510 inline typename tools::promote_args<RT1, RT2, RT3>::type
1511 betac(RT1 a, RT2 b, RT3 x)
1513 return boost::math::betac(a, b, x, policies::policy<>());
1516 template <class RT1, class RT2, class RT3, class Policy>
1517 inline typename tools::promote_args<RT1, RT2, RT3>::type
1518 ibeta(RT1 a, RT2 b, RT3 x, const Policy&)
1520 BOOST_FPU_EXCEPTION_GUARD
1521 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1522 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1523 typedef typename policies::normalise<
1525 policies::promote_float<false>,
1526 policies::promote_double<false>,
1527 policies::discrete_quantile<>,
1528 policies::assert_undefined<> >::type forwarding_policy;
1530 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%)");
1532 template <class RT1, class RT2, class RT3>
1533 inline typename tools::promote_args<RT1, RT2, RT3>::type
1534 ibeta(RT1 a, RT2 b, RT3 x)
1536 return boost::math::ibeta(a, b, x, policies::policy<>());
1539 template <class RT1, class RT2, class RT3, class Policy>
1540 inline typename tools::promote_args<RT1, RT2, RT3>::type
1541 ibetac(RT1 a, RT2 b, RT3 x, const Policy&)
1543 BOOST_FPU_EXCEPTION_GUARD
1544 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1545 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1546 typedef typename policies::normalise<
1548 policies::promote_float<false>,
1549 policies::promote_double<false>,
1550 policies::discrete_quantile<>,
1551 policies::assert_undefined<> >::type forwarding_policy;
1553 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%)");
1555 template <class RT1, class RT2, class RT3>
1556 inline typename tools::promote_args<RT1, RT2, RT3>::type
1557 ibetac(RT1 a, RT2 b, RT3 x)
1559 return boost::math::ibetac(a, b, x, policies::policy<>());
1562 template <class RT1, class RT2, class RT3, class Policy>
1563 inline typename tools::promote_args<RT1, RT2, RT3>::type
1564 ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy&)
1566 BOOST_FPU_EXCEPTION_GUARD
1567 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1568 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1569 typedef typename policies::normalise<
1571 policies::promote_float<false>,
1572 policies::promote_double<false>,
1573 policies::discrete_quantile<>,
1574 policies::assert_undefined<> >::type forwarding_policy;
1576 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%)");
1578 template <class RT1, class RT2, class RT3>
1579 inline typename tools::promote_args<RT1, RT2, RT3>::type
1580 ibeta_derivative(RT1 a, RT2 b, RT3 x)
1582 return boost::math::ibeta_derivative(a, b, x, policies::policy<>());
1586 } // namespace boost
1588 #include <boost/math/special_functions/detail/ibeta_inverse.hpp>
1589 #include <boost/math/special_functions/detail/ibeta_inv_ab.hpp>
1591 #endif // BOOST_MATH_SPECIAL_BETA_HPP