summaryrefslogtreecommitdiff log msg author committer range
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
Diffstat (limited to 'boost/math/special_functions/beta.hpp')
-rw-r--r--boost/math/special_functions/beta.hpp57
1 files changed, 48 insertions, 9 deletions
 diff --git a/boost/math/special_functions/beta.hpp b/boost/math/special_functions/beta.hppindex 8486b82..98d8f7f 100644--- a/boost/math/special_functions/beta.hpp+++ b/boost/math/special_functions/beta.hpp@@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -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 -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()) {- 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(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(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(0)); fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast(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) {