summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/beta.hpp
diff options
context:
space:
mode:
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.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)
{