diff options
Diffstat (limited to 'boost/math/special_functions/beta.hpp')
-rw-r--r-- | boost/math/special_functions/beta.hpp | 57 |
1 files changed, 48 insertions, 9 deletions
diff --git a/boost/math/special_functions/beta.hpp b/boost/math/special_functions/beta.hpp index 8486b824f2..98d8f7fa80 100644 --- a/boost/math/special_functions/beta.hpp +++ b/boost/math/special_functions/beta.hpp @@ -13,6 +13,7 @@ #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/tools/config.hpp> #include <boost/math/special_functions/gamma.hpp> +#include <boost/math/special_functions/binomial.hpp> #include <boost/math/special_functions/factorials.hpp> #include <boost/math/special_functions/erf.hpp> #include <boost/math/special_functions/log1p.hpp> @@ -845,15 +846,54 @@ T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& po // complement of the binomial distribution cdf and use this finite sum. // template <class T> -inline T binomial_ccdf(T n, T k, T x, T y) +T binomial_ccdf(T n, T k, T x, T y) { BOOST_MATH_STD_USING // ADL of std names + T result = pow(x, n); - T term = result; - for(unsigned i = itrunc(T(n - 1)); i > k; --i) + + if(result > tools::min_value<T>()) { - term *= ((i + 1) * y) / ((n - i) * x) ; - result += term; + T term = result; + for(unsigned i = itrunc(T(n - 1)); i > k; --i) + { + term *= ((i + 1) * y) / ((n - i) * x); + result += term; + } + } + else + { + // First term underflows so we need to start at the mode of the + // distribution and work outwards: + int start = itrunc(n * x); + if(start <= k + 1) + start = itrunc(k + 2); + result = pow(x, start) * pow(y, n - start) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start)); + if(result == 0) + { + // OK, starting slightly above the mode didn't work, + // we'll have to sum the terms the old fashioned way: + for(unsigned i = start - 1; i > k; --i) + { + result += pow(x, (int)i) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i)); + } + } + else + { + T term = result; + T start_term = result; + for(unsigned i = start - 1; i > k; --i) + { + term *= ((i + 1) * y) / ((n - i) * x); + result += term; + } + term = start_term; + for(unsigned i = start + 1; i <= n; ++i) + { + term *= (n - i + 1) * x / (i * y); + result += term; + } + } } return result; @@ -1208,9 +1248,9 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de } else if(normalised) { - // the formula here for the non-normalised case is tricky to figure + // The formula here for the non-normalised case is tricky to figure // out (for me!!), and requires two pochhammer calculations rather - // than one, so leave it for now.... + // than one, so leave it for now and only use this in the normalized case.... int n = itrunc(T(floor(b)), pol); T bbar = b - n; if(bbar <= 0) @@ -1221,8 +1261,7 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0)); fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(0)); if(invert) - fract -= (normalised ? 1 : boost::math::beta(a, b, pol)); - //fract = ibeta_series(a+20, bbar, x, fract, l, normalised, p_derivative, y); + fract -= 1; // Note this line would need changing if we ever enable this branch in non-normalized case fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised); if(invert) { |