Probability distribution for fraction of X, given the total sample size N, and the number of Xs in the sample, n. Follows the Appendix of Burgasser et al. (ApJ, 586:512, 2003)

binom_fraction.py

import numpy as np
from scipy.optimize import bisect
from scipy.special import binom as binom_coeff
from scipy.integrate import quad
from functools import partial

def binom_function(N, n, p):
	c = binom_coeff(N, n)
	return c * p**n * (1.-p)**(N-n)

def Bprime(N, n, p):
	""" This is B' from Burgasser+(2003) """
	B = binom_function(N, n, p)
	return (N+1)*B

def integrate_upto(N, n, up):
	f = partial(Bprime, N, n)
	return quad(f, 0, up)

def integrate_from(N, n, low):
	f = partial(Bprime, N, n)
	return quad(f, low, 1.)


def low_upp_limits(n, N):
	if n==0:
		ll = 0.
		ul = bisect(lambda p: integrate_upto(N, n, p)[0] - 0.84, 0, 1) * 100  
	else:
		ll = bisect(lambda p: integrate_from(N, n, p)[0] - 0.84, 0, 1) * 100
		ul = bisect(lambda p: integrate_upto(N, n, p)[0] - 0.84, 0, 1) * 100
    
	print 'lower limit', '%3.2f'% (ll), '%'
	print 'upper limit', '%3.2f'% (ul), '%'

	return ll, ul

def calc_mode(n, N):
  """ The mode of the B' distribution (same as for the B distribution) """
	return (1+n-1)/(2.+N-2)*100



if __name__ == '__main__':
	# example from Burgasser+(2003)
  N = 10
  n = 2

  mode = calc_mode(n, N)
  print 'mode', '%3.2f' % (mode), '%'
  print '\n\n'

  ll, ul = low_upp_limits(n, N)
  print '$%3.2f\%%^{+%3.2f\%%}_{-%3.2f\%%}$' % (mode, ul-mode, mode-ll)