# Computational shortcuts/approximations for beta-binomial and beta-negative binomial CDFs

There are several simplifications that can be done so that computing cumulative distribution functions of beta-binomial and beta-negative binomial distributions, but still computing CDF as $F(x) = \sum_{k=0}^x f(k)$ is computationally intensive. Are there any computational shortcuts, or approximations that can be used so to simplify the calculations?

• It seems you miss something, should it be$F(x)= \sum_{k=0}^x f(k)\Delta k$? – Deep North May 16 '17 at 1:02
• @DeepNorth - the two distributions are discrete, more precisely, on the non-negative integers. – jbowman May 16 '17 at 3:59
• I have some code at iandjmsmith.wordpress.com for both. It's for Excel spreadsheets but I also have other code in Pascal for both, if you would prefer. Ian Smith – Ian Smith May 30 '17 at 12:44

I have no idea what you know on this subject which is not particularly simple.

Obviously, you can express f(k+1) in terms of f(k) to save the effort of calculating the pdf function x+1 times. To minimise error propagation, you would be best finding the largest value of f(k) and working with summations away from that value until the terms are sufficiently small. If x is large, say > 100million then this approach may still not be very good. You also cannot use this approach and subtract from 1 to find the summation from x+1 upwards because the summation to x can be arbitrarily close to 1.

On the other hand if you only want to calculate the cdf and only for small values of x, then this type of approach will be adequate.

Ian Smith

P.S. There are obvious connections between cdf_hypergeometric, cdf_betabinomial & cdf_betanegativebinomial when the parameter values are all integral but they can be extended when the parameters are non-integral.

cdf_BNB(i,r,a,b) = pmf_BNB(i,r,b,a) + cdf_BNB(i-1,r,a,b)

cdf_BNB(i,r,a,b) = PBB(i,r,b,a) + cdf_BNB(i-1,r+1,a,b)
and since cdf_BNB(i,r,a,b) = cdf_BNB(i,b,a,r)
cdf_BNB(i,r,a,b) = PBB (i,b,r,a) + cdf_BNB(i-1,r,a,b+1)

cdf_BNB(i,r,a,b) = HT(i,r,b,a-1) + cdf_BNB(i-1,r+1,a-1,b+1)

where PBB is essentially the pmf of the BetaBinomial but where non-integral values are allowed for the sample size i.e. PBB(i,r,a,b) = pmf_BB(i,i+r,a,b) when r is integral and HT is essentially the pmf of the hypergeometric distribution but where non-integral values are allowed for the sample size, total number of type1s and total number of items to be selected from. HT(i,r,a,b) = pmf_hypergeometric(i,i+r,i+a,i+r+a+b)

The equivalent survival function relationships are

sf_BNB(i,r,a,b) = pmf_BNB(i+1,r,a,b) + sf_BNB(i+1,r,a,b)

sf_BNB(i,r,a,b) = PBB(i+1,r-1,a,b) + sf_BNB(i+1,r-1,a,b)
sf_BNB(i,r,a,b) = PBB(i+1,b-1,a,r) + sf_BNB(i+1,r,a,b-1)
and

sf_BNB(i,r,a,b) = HT(i+1,r-1,b-1,a) + sf_BNB(i+1,r-1,a+1,b-1)

For integral r, b
sf_BNB(i,r,a,b) =cdf_BNB(r-1,i+1,b,a)
sf_BNB(i,r,a,b) =cdf_BNB(b-1,i+1,r,a)

The python code in mystatsBetaBinomialpy.txt in https://github.com/scipy/scipy/issues/7102 will calculate the cdf and sf (cumulative distribution function & survival function) of the Beta Binomial. It also calculates the cdfs and sfs of the Beta Negative Binomial and the hypergeometric distribution which are all closely related. The python code in mystatsFetpy.txt in https://github.com/scipy/scipy/issues/9231 continues with another use of the identities used in the algorithms for the cdf & sf of the Beta Negative Binomial distribution. The code is for Fisher’s exact test for 2x2, 2x3 & 2x4 tables. These are designed to work quickly with tables with large entries.

The following are some brief notes on how the code for the Beta Binomial works. The cdf and sf of the Beta Binomial can be expressed in terms of the sf of the Beta Negative Binomial. Basically,

cdf_BB(i, sample_size, beta_shape1, beta_shape2) = sf_BNB (sample_size – i - 1, i + 1, beta_shape1, beta_shape2)


and

sf_BB(i, sample_size, beta_shape1, beta_shape2) = 0 if i >= sample_size, otherwise

sf_BB(i, sample_size, beta_shape1, beta_shape2) = sf_BNB (i, sample_size - i, beta_shape2, beta_shape2).


The main driver of accuracy in cdf and sf calculations is usually the accuracy of the pmf functions. Often implementations of the pmf function involving factorials call functions delivering the log of the gamma function, sum a few of them and exponentiate the result. Unfortunately, gammaln(large number) returns a larger number and

exp(sum of large numbers)

is exp(answer with large absolute error)

which returns answer with large relative error.

This can be virtually cured by using functions which deliver the log of the gamma function in quadruple precision where “sum of large numbers” will deliver “answer with small absolute error” or mitigated against by using functions which deliver the log of the gamma function in long double precision (80-bit arithmetic) but this is only a case of getting 3-4 more significant figures.

Assuming we can’t or don’t want to do this, the python code I have referred to, does pretty much as well as you can with double precision, i.e. the relative error is roughly ln(true pmf value).1e-16. The code uses the pmf function for the hypergeometric and details of how that function works are given at the end of these notes.

In the other reply I gave, I listed a few of the properties which can be used to perform the calculations of the cdf & sf functions. If you look at the main else section of the sf_BetaNegativeBinomial code, it uses all of the following sf relationships.

sf_BNB(i,r,a,b) = pmf_BNB(i+1,r,a,b) + sf_BNB(i+1,r,a,b)

sf_BNB(i,r,a,b) = PBB(i+1,r-1,b,a) + sf_BNB(i+1,r-1,a,b)

sf_BNB(i,r,a,b) = PBB(i+1,b-1,r,a) + sf_BNB(i+1,r,a,b-1)

sf_BNB(i,r,a,b) = HT(i+1,r-1,b-1,a) + sf_BNB(i+1,r-1,a+1,b-1)

sf_BNB(i,r,a,b) = PBB(a,r-1,b, i+1) + sf_BNB(i,r-1,a+1,b)

sf_BNB(i,r,a,b) = PBB(a,b-1,r, i+1) + sf_BNB(i,r,a+1,b-1)

sf_BNB(i,r,a,b) = pBNB(a,r,i+1,b) + sf_BNB(i,r,a+1,b)


where pBNB is the same as pmf_BNB except the check for the first parameter being integral has been removed, PBB is the same as pmf_BB except that checks have been removed and pmf_BB(i,ss,a,b) = PBB(i,ss-i,a,b) and pmf_hypergeometric(t1,ss,total_t1,ps) = HT(t1,ss-t1,total_t1-t1,ps-total_t1-ss+t1). Again checks are removed on the parameters of HT so the function can be easily used to evaluate pmf_BNB.

All this still isn’t adequate to deal with the full range of the parameters. There are always problems with summation algorithms, particularly with large parameter values. The continued fraction algorithm for the hypergeometric distribution is the most important technique for taking care of large parameter values. Details of how the hypergeometric distribution functions work are given at the end of these notes. That still leaves one area where calculations take a long time because large numbers of summations have to be done before the continued fraction algorithm can be used. The examples sf_BetaNegativeBinomial(1e6,1e14,1e-5,1e-14) and cdf_BetaNegativeBinomial(4503599627370495,30000000.5,0.005,30000000.5) in the Beta Binomial code are included to show how slow things can be and give ideas for choosing limits on i+sample_size+beta_shape1+beta_shape2.

As the BetaBinomial, BetaNegativeBinomial and Hypergeometric distributions are closely connected, the sf relationships above can be transformed for application to the sf and cdf functions for the Hypergeometric distribution. An example which is easy to verify is

cdf_hypergeometric (type1s, sample_size, total_type1s, pop_size)
+ pmf_BB(pop_size - total_type1s - sample_size + type1s +1, sample_size - type1s -1, total_type1s - type1s,  type1s + 1)
= cdf_hypergeometric (type1s, sample_size-1, total_type1s, pop_size)


Such relationships involving the cdf & sf functions for the Hypergeometric distribution are not of any great value I know of, except when it comes to Fisher’s exact test. With the 2x2 Fisher’s exact test, the probability of a given table with entries a,b,c,d is pmf_hypergeometric(a,a+b,a+c,a+b+c+d). If we consider the set of tables a+e,b-e,c-e,d+e for all integer values of e then these are all the tables with fixed row totals a+b,c+d and fixed column totals a+c, b+d. So summing the probabilities of all tables whose probabilities are as small or smaller than the observed table (a,b,c,d) will give us the Fisher’s exact test value.

The probability of each table is pmf_hypergeometric(a+e, a+b, a+c, a+b+c+d) and so the value for a 2-sided Fisher’s exact test can be expressed as cdf_hypergeometric(a+e0, a+b, a+c, a+b+c+d) + sf_hypergeometric(a+e1, a+b, a+c, a+b+c+d) where one or both of the values e0, e1 will be 0.

For a 2x3 Fisher’s test the probability of the table with rows [a,b,c],[d,e,f] can be calculated as pmf_hypergeometric(a, a+b, a+d, a+b+d+e). pmf_hypergeometric(c, c+f, a+b+c, a+b+c+d+e+f) and if we consider the set of all tables with rows of the form [a+g,b-g+h,c-h],[d-g,e+g-h,f+h] where g & h are integers then these are all the tables with fixed row totals a+b+c, d+e+f and fixed column totals a+d, b+e, c+f. The table probabilities can be calculated as pmf_hypergeometric(a+g, a+b+h, a+d, a+b+d+e). pmf_hypergeometric(c-h, c+f, a+b+c, a+b+c+d+e+f) So if we fix h first and sum over all possible values of g the total is (cdf_hypergeometric(a+g0(h), a+b+h, a+d, a+b+d+e) + sf_hypergeometric(a+g1(h), a+b+h, a+d, a+b+d+e)).pmf_hypergeometric(c-h, c+f, a+b+c, a+b+c+d+e+f) and summing over the possible values of h gives us the answer.

The main shortcut is that we can express the cdf_hypergeometric(a+g0(h), a+b+h-1, a+d, a+b+d+e) as cdf_hypergeometric(a+g0(h), a+b+h, a+d, a+b+d+e) plus a pmf_function, so potentially time-consuming evaluations of the cdf and sf functions are avoided. Of course, we still have to find out what g0(h-1) and g1(h-1) are and adjust the cdf and sf values for this but it is a lot quicker than finding g0(h-1) and g1(h-1) and the calculating the cdf and sf functions from scratch. There are also another couple of shortcuts which can be applied. First if we choose the initial value of h to be one where a mode is known to occur and the values of g0(hmode) and g1(hmode) mean all values are to be included then the table with the largest probability is to be included and so all tables must be included and so the answer is 1. If we work our way out from the value hmode and some g0(h), g1(h) mean that all tables are to be included for that value of h, then all tables for more extreme values of h must also be included.

For a 2x4 Fisher’s test the probability of the table with rows [a,b,c,d],[e,f,g,h] can be calculated as pmf_hypergeometric(a, a+b, a+e, a+b+e+f).pmf_hypergeometric(a+b, a+b+e+f, a+b+c+d, a+b+c+d+e+f+g+h).pmf_hypergeometric(d, d+h, c+d, c+d+g+h) and if we consider the set of all tables with rows of the form [a+i,b-i+j,c-j+k,d-k], [e-i,f+i-j,g+j-k,h+k] where i, j & k are integers then these are all the tables with fixed row totals a+b+c+d, e+f+g+h and fixed column totals a+e, b+f, c+g, d+h. The table probabilities can be calculated as pmf_hypergeometric(a+i, a+b+j, a+e, a+b+e+f) . pmf_hypergeometric(a+b+j, a+b+e+f, a+b+c+d, a+b+c+d+e+f+g+h) . pmf_hypergeometric(d-k, d+h, c+d-j, c+d+g+h) Here, we fix j first, then for all k we can find suitable values of i so we can sum (cdf_hypergeometric(a+i0(k,j), a+b+j, a+e, a+b+e+f) + sf_hypergeometric(a+i1(k,j), a+b+j, a+e, a+b+e+f) ) .pmf_hypergeometric(d-k, d+h, c+d-j, c+d+g+h) over k. Finally, we sum these over all j. Note that if we arrange to run through all possible k in decreasing order of pmf_hypergeometric(d-k, d+h, c+d-j, c+d+g+h), then a+i0(k,j) is non-decreasing and a+i1(k,j) is non-increasing with the decreasing pmf values. So the cdf and sf functions are easy to modify with k and we can use our rules to change the sample size parameter for a hypergeometric distribution to recalculate the cdf and sf functions when we change j. Running through all possible k in decreasing order of pmf_hypergeometric(d-k, d+h, c+d-j, c+d+g+h) is trivial. We start at the modal value of k and the next biggest pmf value is either at k+1 or k -1 and we continue working our way out.

We can use the same factorization scheme for 2x5 tables, unfortunately, running through all possible 2x3 tables in decreasing order of pmf becomes a problem. I have methods for evaluating Fisher’s exact for 2x5 tables but I cannot say these are the best possible for tables with large entries.

Notes on Hypergeometric distribution functions.

If we have a population size of m items, where k, k <= m, of the objects have a certain attribute and all m items have an equal chance of being selected and we select n items, n>= 0 & n <= m, without replacement, then the random variable X, which is the number of items selected with the attribute, has a hypergeometric distribution.

Pr(X=i;n,k,m)
= (m-n)!n!(m-k)!k!/i!/(n-i)!/(k-i)!/(m-k-n+i)!/m! if i >= 0 and i <= n and i <= k and k+n <= m+i
= 0 otherwise.


Using i! = sqrt(2.pi).exp(-(i+1)+logfbit(i)).(i+1)**(i+½)

Pr(X=i;n,k,m)

= exp(1+logfbit(m-n)+logfbit(n)+logfbit(m-k)+logfbit(k)-logfbit(i)-logfbit(n-i)-logfbit(k-i)-logfbit(m-k-n+i)-logfbit(m))

.((m-n+1)(m-k+1)/(m-k-n+i+1)/(m+1))**(m-n-k+i)

.((k+1)(m-n+1)/(k-i+1)/(m+1))**(k-i)

.((n+1)(m-k+1)/(n-i+1)/(m+1))**(n-i)

.((k+1)(n+1)/(i+1)/(m+1))i

/sqrt(2.pi.(k+1)(m-k+1)(n+1)(m-n+1)/(m-k-n+i+1)/(k-i+1)/(n-i+1)/(i+1)/(m+1))

if i >= 0 and i <= n and i <= k and k+n <= m+i

= 1 if i = 0

= 0 otherwise.


Defining log1(x) to be ln(1+x)-x

we have ln(((k+1).(n+1)/(i+1)/(m+1))i)

== i.ln(1+ ((k+1).(n+1)/(i+1)/(m+1)-1)

== i.(ln(1+(m.i-n.k-i)/(i+1)/(m+1)) - (m.i-n.k-i)/(i+1)/(m+1)) - (m.i-n.k-i)/(i+1)/(m+1) + (m.i-n.k-i)/(m+1)

== i.log1((m.i-n.k-i)/(i+1)/(m+1)) - (m.i-n.k-i)/(i+1)/(m+1) + (m.i-n.k-i)/(m+1)


With similar simplifications for the other powers, we get

 ln(

((m-n+1)(m-k+1)/(m-k-n+i+1)/(m+1))(m-n-k+i)

.((k+1)(m-n+1)/(k-i+1)/(m+1))(k-i)

.((n+1)(m-k+1)/(n-i+1)/(m+1))(n-i)

.((k+1)(n+1)/(i+1)/(m+1))i

)


which is

(m-n-k+i).log1((m.i-n.k-i)/(m-k-n+i+1)/(m+1)) - (m.i-n.k-i)/(i+1)/(m+1) + (m.i-n.k-i)/(m+1)

+(k-i).log1((-(m.i-n.k)-(n-i))/(k-i+1)/(m+1)) - (m.i-n.k-(n-i))/(k-i+1)/(m+1) + (-(m.i-n.k)-(n-i))/(m+1)

+(n-i).log1((-(m.i-n.k)-(k-i))/(n-i+1)/(m+1)) - (m.i-n.k-(k-i))/(n-i+1)/(m+1)+ (-(m.i-n.k)-(k-i))/(m+1)

+i.log1((m.i-n.k-(m-k-n+i))/(i+1)/(m+1)) - (m.i-n.k-(m-k-n+i))/(i+1)/(m+1) + (m.i-n.k-(m-k-n+i))/(m+1)


which is

(m-n-k+i).log1((m.i-n.k-i)/(m-k-n+i+1)/(m+1)) - (m.i-n.k-i)/(i+1)/(m+1)

+(k-i).log1((-(m.i-n.k)-(n-i))/(k-i+1)/(m+1)) - (m.i-n.k-(n-i))/(k-i +1)/(m+1)

+(n-i).log1((-(m.i-n.k)-(k-i))/(n-i+1)/(m+1)) - (m.i-n.k-(k-i))/(n-i+1)/(m+1)

+i.log1((m.i-n.k-(m-k-n+i))/(i+1)/(m+1)) - (m.i-n.k-(m-k-n+i))/(i+1)/(m+1)
-1+1/m


Giving us

Pr(X=i;n,k,m)

= exp(
1/m+logfbit(m-n)+logfbit(n)+logfbit(m-k)+logfbit(k)-logfbit(i)-logfbit(n-i)-logfbit(k-i)-logfbit(m-k-n+i)-logfbit(m)

+(m-n-k+i).log1((m.i-n.k-i)/(m-k-n+i+1)/(m+1)) - (m.i-n.k-i)/(i+1)/(m+1)

+(k-i).log1((-(m.i-n.k)-(n-i))/(k-i+1)/(m+1)) - (m.i-n.k-(n-i))/(k-i +1)/(m+1)

+(n-i).log1((-(m.i-n.k)-(k-i))/(n-i+1)/(m+1)) - (m.i-n.k-(k-i))/(n-i+1)/(m+1)

+i.log1((m.i-n.k-(m-k-n+i))/(i+1)/(m+1)) - (m.i-n.k-(m-k-n+i))/(i+1)/(m+1)

)/sqrt(2.pi.(k+1)(m-k+1)(n+1)(m-n+1)/(m-k-n+i+1)/(k-i+1)/(n-i+1)/(i+1)/(m+1))


Note for large values for i,k,n & m, the most important term for accuracy in the pmf calculation is m.i-n.k The code uses a function Generalabminuscd which evaluates this expression accurately. My lack of experience with python means that this is a function which might be easily improved. Note also that if the arguments of log1 are too near -1, we use the original expressions such as i.ln((k+1)(n+1)/(i+1)/(m+1)) instead of i.log1((m.i-n.k-(m-k-n+i))/(i+1)/(m+1))

The cumulative distribution function for the hypergeometric distribution is evaluated using the following equations

For (i+½).(m+½) <= (n-½).(k-½)

Pr(X<=i;n,k,m) = cf(i).Pr(X=i;n,k,m)


where cf(i) = 1 + cf(i-1).i.(m-k-n+i) /(n-i+1)/(k-i+1) {Summation Equation }

cf(i)

= 1 + i.(m-k-n+i)/

((n-i+1).(k-i+1)-(i-1).(m-k-n+i-1) + m.(i-1).(m-k-n+i-1)/

((n-i+2).(k-i+2)-(i-3).(m-k-n+i-2)+(i-1) + 2.(m-1).(i-2).(m-k-n+i-2)/

((n-i+3).(k-i+3)-(i-5).(m-k-n+i-3)+2.(i-2) + 3.(m-2).(i-3).(m-k-n+i-3)/

((n-i+4).(k-i+4)-(i-7).(m-k-n+i-4)+3.(i-3) + ...) {Continued Fraction Equation}


When i is reasonably large and close to (n-½).(k-½)/(m+½) the continued fraction equation can be used to sum the tail of the distribution (i.e. we use it to evaluate cf(i-x) for some x) and then the summation equation used to find cf(i).

For (i+½).(m+½) > (n-½).(k-½)

Pr(X<=i;n,k,m) = 1-Pr(X<=(n-i-1);n,m-k,m)

Ian Smith

• Please, this needs some better explanations, for instance, what is the sf_ functions? Also needs better markup. I started a cleanup, but please complete, at least by marking all the code. – kjetil b halvorsen May 19 at 14:02