<s> """<STR_LIT>""" <EOL> from math import exp , log , sqrt <EOL> LSTEP = <NUM_LIT:200> <EOL> EXPTHRES = exp ( LSTEP ) <EOL> EXPSTEP = exp ( - LSTEP ) <EOL> def normal_01_cdf ( x ) : <EOL> """<STR_LIT>""" <EOL> a1 = <NUM_LIT> <EOL> a2 = <NUM_LIT> <EOL> a3 = <NUM_LIT> <EOL> a4 = <NUM_LIT> <EOL> a5 = <NUM_LIT> <EOL> a6 = <NUM_LIT> <EOL> a7 = <NUM_LIT> <EOL> b0 = <NUM_LIT> <EOL> b1 = <NUM_LIT> <EOL> b2 = <NUM_LIT> <EOL> b3 = <NUM_LIT> <EOL> b4 = <NUM_LIT> <EOL> b5 = <NUM_LIT> <EOL> b6 = <NUM_LIT> <EOL> b7 = <NUM_LIT> <EOL> b8 = <NUM_LIT> <EOL> b9 = <NUM_LIT> <EOL> b10 = <NUM_LIT> <EOL> b11 = <NUM_LIT> <EOL> if abs ( x ) <= <NUM_LIT> : <EOL> y = <NUM_LIT:0.5> * x * x <EOL> q = <NUM_LIT:0.5> - abs ( x ) * ( a1 - a2 * y / ( y + a3 - a4 / ( y + a5 + a6 / ( y + a7 ) ) ) ) <EOL> elif abs ( x ) <= <NUM_LIT> : <EOL> y = <NUM_LIT:0.5> * x * x <EOL> q = exp ( - y ) * b0 / ( abs ( x ) - b1 + b2 / ( abs ( x ) + b3 + b4 / ( abs ( x ) - b5 + b6 / ( abs ( x ) + b7 - b8 / ( abs ( x ) + b9 + b10 / ( abs ( x ) + b11 ) ) ) ) ) ) <EOL> else : <EOL> q = <NUM_LIT:0.0> <EOL> if x < <NUM_LIT:0.0> : <EOL> cdf = q <EOL> else : <EOL> cdf = <NUM_LIT:1.0> - q <EOL> return cdf <EOL> def normal_cdf_inv ( p , mu = None , sigma2 = None , lower = True ) : <EOL> upper = not lower <EOL> if p < <NUM_LIT:0> or p > <NUM_LIT:1> : <EOL> raise Exception ( "<STR_LIT>" % p ) <EOL> split = <NUM_LIT> <EOL> a0 = <NUM_LIT> <EOL> a1 = - <NUM_LIT> <EOL> a2 = <NUM_LIT> <EOL> a3 = - <NUM_LIT> <EOL> b1 = - <NUM_LIT> <EOL> b2 = <NUM_LIT> <EOL> b3 = - <NUM_LIT> <EOL> b4 = <NUM_LIT> <EOL> c0 = - <NUM_LIT> <EOL> c1 = - <NUM_LIT> <EOL> c2 = <NUM_LIT> <EOL> c3 = <NUM_LIT> <EOL> d1 = <NUM_LIT> <EOL> d2 = <NUM_LIT> <EOL> q = p - <NUM_LIT:0.5> <EOL> r = <NUM_LIT:0.0> <EOL> ppnd = <NUM_LIT:0.0> <EOL> if abs ( q ) <= split : <EOL> r = q * q <EOL> ppnd = q * ( ( ( a3 * r + a2 ) * r + a1 ) * r + a0 ) / ( ( ( ( b4 * r + b3 ) * r + b2 ) * r + b1 ) * r + <NUM_LIT:1> ) <EOL> else : <EOL> r = p <EOL> if q > <NUM_LIT:0> : <EOL> r = <NUM_LIT:1> - p <EOL> if r > <NUM_LIT:0> : <EOL> r = sqrt ( - log ( r ) ) <EOL> ppnd = ( ( ( c3 * r + c2 ) * r + c1 ) * r + c0 ) / ( ( d2 * r + d1 ) * r + <NUM_LIT:1> ) <EOL> if q < <NUM_LIT:0> : <EOL> ppnd = - ppnd <EOL> else : <EOL> ppnd = <NUM_LIT:0> <EOL> if upper : <EOL> ppnd = - ppnd <EOL> if mu != None and sigma2 != None : <EOL> return ppnd * sqrt ( sigma2 ) + mu <EOL> else : <EOL> return ppnd <EOL> def normal_cdf ( z , mu = <NUM_LIT:0.0> , sigma2 = <NUM_LIT:1.0> , lower = True ) : <EOL> upper = not lower <EOL> z = ( z - mu ) / sqrt ( sigma2 ) <EOL> ltone = <NUM_LIT> <EOL> utzero = <NUM_LIT> <EOL> con = <NUM_LIT> <EOL> a1 = <NUM_LIT> <EOL> a2 = <NUM_LIT> <EOL> a3 = <NUM_LIT> <EOL> a4 = <NUM_LIT> <EOL> a5 = <NUM_LIT> <EOL> a6 = <NUM_LIT> <EOL> a7 = <NUM_LIT> <EOL> b1 = <NUM_LIT> <EOL> b2 = <NUM_LIT> <EOL> b3 = <NUM_LIT> <EOL> b4 = <NUM_LIT> <EOL> b5 = <NUM_LIT> <EOL> b6 = <NUM_LIT> <EOL> b7 = <NUM_LIT> <EOL> b8 = <NUM_LIT> <EOL> b9 = <NUM_LIT> <EOL> b10 = <NUM_LIT> <EOL> b11 = <NUM_LIT> <EOL> b12 = <NUM_LIT> <EOL> y = <NUM_LIT:0.0> <EOL> alnorm = <NUM_LIT:0.0> <EOL> if z < <NUM_LIT:0> : <EOL> upper = not upper <EOL> z = - z <EOL> if z <= ltone or upper and z <= utzero : <EOL> y = <NUM_LIT:0.5> * z * z <EOL> if z > con : <EOL> alnorm = b1 * exp ( - y ) / ( z - b2 + b3 / ( z + b4 + b5 / ( z - b6 + b7 / ( z + b8 - b9 / ( z + b10 + b11 / ( z + b12 ) ) ) ) ) ) <EOL> else : <EOL> alnorm = <NUM_LIT:0.5> - z * ( a1 - a2 * y / ( y + a3 - a4 / ( y + a5 + a6 / ( y + a7 ) ) ) ) <EOL> else : <EOL> alnorm = <NUM_LIT:0.0> <EOL> if not upper : <EOL> alnorm = <NUM_LIT:1.0> - alnorm <EOL> return alnorm <EOL> def poisson_cdf ( n , lam , lower = True ) : <EOL> """<STR_LIT>""" <EOL> k = int ( n ) <EOL> if lam <= <NUM_LIT:0.0> : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> if lower : <EOL> if lam > <NUM_LIT> : <EOL> return __poisson_cdf_large_lambda ( k , lam ) <EOL> else : <EOL> return __poisson_cdf ( k , lam ) <EOL> else : <EOL> if lam > <NUM_LIT> : <EOL> return __poisson_cdf_Q_large_lambda ( k , lam ) <EOL> else : <EOL> return __poisson_cdf_Q ( k , lam ) <EOL> def __poisson_cdf ( k , a ) : <EOL> """<STR_LIT>""" <EOL> if k < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> next = exp ( - a ) <EOL> cdf = next <EOL> for i in xrange ( <NUM_LIT:1> , k + <NUM_LIT:1> ) : <EOL> last = next <EOL> next = last * a / i <EOL> cdf = cdf + next <EOL> if cdf > <NUM_LIT:1> : <EOL> return <NUM_LIT:1> <EOL> else : <EOL> return cdf <EOL> def __poisson_cdf_large_lambda ( k , a ) : <EOL> """<STR_LIT>""" <EOL> if k < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> num_parts = int ( a / LSTEP ) <EOL> last_part = a % LSTEP <EOL> lastexp = exp ( - last_part ) <EOL> next = EXPSTEP <EOL> num_parts -= <NUM_LIT:1> <EOL> cdf = next <EOL> for i in xrange ( <NUM_LIT:1> , k + <NUM_LIT:1> ) : <EOL> last = next <EOL> next = last * a / i <EOL> cdf = cdf + next <EOL> if next > EXPTHRES or cdf > EXPTHRES : <EOL> if num_parts >= <NUM_LIT:1> : <EOL> cdf *= EXPSTEP <EOL> next *= EXPSTEP <EOL> num_parts -= <NUM_LIT:1> <EOL> else : <EOL> cdf *= lastexp <EOL> lastexp = <NUM_LIT:1> <EOL> for i in xrange ( num_parts ) : <EOL> cdf *= EXPSTEP <EOL> cdf *= lastexp <EOL> return cdf <EOL> def __poisson_cdf_Q ( k , a ) : <EOL> """<STR_LIT>""" <EOL> if k < <NUM_LIT:0> : <EOL> return <NUM_LIT:1> <EOL> next = exp ( - a ) <EOL> for i in xrange ( <NUM_LIT:1> , k + <NUM_LIT:1> ) : <EOL> last = next <EOL> next = last * a / i <EOL> cdf = <NUM_LIT:0> <EOL> i = k + <NUM_LIT:1> <EOL> while next > <NUM_LIT:0> : <EOL> last = next <EOL> next = last * a / i <EOL> cdf += next <EOL> i += <NUM_LIT:1> <EOL> return cdf <EOL> def __poisson_cdf_Q_large_lambda ( k , a ) : <EOL> """<STR_LIT>""" <EOL> if k < <NUM_LIT:0> : <EOL> return <NUM_LIT:1> <EOL> num_parts = int ( a / LSTEP ) <EOL> last_part = a % LSTEP <EOL> lastexp = exp ( - last_part ) <EOL> next = EXPSTEP <EOL> num_parts -= <NUM_LIT:1> <EOL> for i in xrange ( <NUM_LIT:1> , k + <NUM_LIT:1> ) : <EOL> last = next <EOL> next = last * a / i <EOL> if next > EXPTHRES : <EOL> if num_parts >= <NUM_LIT:1> : <EOL> next *= EXPSTEP <EOL> num_parts -= <NUM_LIT:1> <EOL> else : <EOL> cdf *= lastexp <EOL> lastexp = <NUM_LIT:1> <EOL> cdf = <NUM_LIT:0> <EOL> i = k + <NUM_LIT:1> <EOL> while next > <NUM_LIT:0> : <EOL> last = next <EOL> next = last * a / i <EOL> cdf += next <EOL> i += <NUM_LIT:1> <EOL> if next > EXPTHRES or cdf > EXPTHRES : <EOL> if num_parts >= <NUM_LIT:1> : <EOL> cdf *= EXPSTEP <EOL> next *= EXPSTEP <EOL> num_parts -= <NUM_LIT:1> <EOL> else : <EOL> cdf *= lastexp <EOL> lastexp = <NUM_LIT:1> <EOL> for i in xrange ( num_parts ) : <EOL> cdf *= EXPSTEP <EOL> cdf *= lastexp <EOL> return cdf <EOL> def poisson_cdf_inv ( cdf , lam , maximum = <NUM_LIT:1000> ) : <EOL> """<STR_LIT>""" <EOL> assert lam < <NUM_LIT> <EOL> if cdf < <NUM_LIT:0> or cdf > <NUM_LIT:1> : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> elif cdf == <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> sum2 = <NUM_LIT:0> <EOL> newval = exp ( - lam ) <EOL> sum2 = newval <EOL> for i in xrange ( <NUM_LIT:1> , maximum + <NUM_LIT:1> ) : <EOL> sumold = sum2 <EOL> last = newval <EOL> newval = last * lam / i <EOL> sum2 = sum2 + newval <EOL> if sumold <= cdf and cdf <= sum2 : <EOL> return i <EOL> return maximum <EOL> def poisson_cdf_Q_inv ( cdf , lam , maximum = <NUM_LIT:1000> ) : <EOL> """<STR_LIT>""" <EOL> assert lam < <NUM_LIT> <EOL> if cdf < <NUM_LIT:0> or cdf > <NUM_LIT:1> : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> elif cdf == <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> sum2 = <NUM_LIT:0> <EOL> newval = exp ( - lam ) <EOL> sum2 = newval <EOL> for i in xrange ( <NUM_LIT:1> , maximum + <NUM_LIT:1> ) : <EOL> sumold = sum2 <EOL> last = newval <EOL> newval = last * lam / i <EOL> sum2 = sum2 + newval <EOL> if sumold <= cdf and cdf <= sum2 : <EOL> return i <EOL> return maximum <EOL> def poisson_pdf ( k , a ) : <EOL> """<STR_LIT>""" <EOL> if a <= <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> return exp ( - a ) * pow ( a , k ) / factorial ( k ) <EOL> def binomial_coef ( n , k ) : <EOL> """<STR_LIT>""" <EOL> mn = min ( k , n - k ) <EOL> if mn < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> elif mn == <NUM_LIT:0> : <EOL> return <NUM_LIT:1> <EOL> else : <EOL> mx = max ( k , n - k ) <EOL> cnk = float ( mx + <NUM_LIT:1> ) <EOL> for i in xrange ( <NUM_LIT:2> , mn + <NUM_LIT:1> ) : <EOL> cnk = cnk * ( mx + i ) / i <EOL> return cnk <EOL> def binomial_cdf ( x , a , b , lower = True ) : <EOL> """<STR_LIT>""" <EOL> if lower : <EOL> return _binomial_cdf_f ( x , a , b ) <EOL> else : <EOL> return _binomial_cdf_r ( x , a , b ) <EOL> def _binomial_cdf_r ( x , a , b ) : <EOL> if x < <NUM_LIT:0> : <EOL> return <NUM_LIT:1> <EOL> elif a < x : <EOL> return <NUM_LIT:0> <EOL> elif b == <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> elif b == <NUM_LIT:1> : <EOL> return <NUM_LIT:1> <EOL> else : <EOL> argmax = int ( a * b ) <EOL> if x < argmax : <EOL> seedpdf = binomial_pdf ( argmax , a , b ) <EOL> pdf = seedpdf <EOL> cdf = pdf <EOL> for i in xrange ( argmax - <NUM_LIT:1> , x , - <NUM_LIT:1> ) : <EOL> pdf /= ( a - i ) * b / ( <NUM_LIT:1> - b ) / ( i + <NUM_LIT:1> ) <EOL> if pdf == <NUM_LIT:0.0> : break <EOL> cdf += pdf <EOL> pdf = seedpdf <EOL> i = argmax <EOL> while True : <EOL> pdf *= ( a - i ) * b / ( <NUM_LIT:1> - b ) / ( i + <NUM_LIT:1> ) <EOL> if pdf == <NUM_LIT:0.0> : break <EOL> cdf += pdf <EOL> i += <NUM_LIT:1> <EOL> cdf = min ( <NUM_LIT:1> , cdf ) <EOL> cdf = float ( "<STR_LIT>" % cdf ) <EOL> return cdf <EOL> else : <EOL> pdf = binomial_pdf ( x + <NUM_LIT:1> , a , b ) <EOL> cdf = pdf <EOL> i = x + <NUM_LIT:1> <EOL> while True : <EOL> pdf *= ( a - i ) * b / ( <NUM_LIT:1> - b ) / ( i + <NUM_LIT:1> ) <EOL> if pdf == <NUM_LIT:0.0> : break <EOL> cdf += pdf <EOL> i += <NUM_LIT:1> <EOL> cdf = min ( <NUM_LIT:1> , cdf ) <EOL> cdf = float ( "<STR_LIT>" % cdf ) <EOL> return cdf <EOL> def _binomial_cdf_f ( x , a , b ) : <EOL> if x < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> elif a < x : <EOL> return <NUM_LIT:1> <EOL> elif b == <NUM_LIT:0> : <EOL> return <NUM_LIT:1> <EOL> elif b == <NUM_LIT:1> : <EOL> return <NUM_LIT:0> <EOL> else : <EOL> argmax = int ( a * b ) <EOL> if x > argmax : <EOL> seedpdf = binomial_pdf ( argmax , a , b ) <EOL> pdf = seedpdf <EOL> cdf = pdf <EOL> for i in xrange ( argmax - <NUM_LIT:1> , - <NUM_LIT:1> , - <NUM_LIT:1> ) : <EOL> pdf /= ( a - i ) * b / ( <NUM_LIT:1> - b ) / ( i + <NUM_LIT:1> ) <EOL> if pdf == <NUM_LIT:0.0> : break <EOL> cdf += pdf <EOL> pdf = seedpdf <EOL> for i in xrange ( argmax , x ) : <EOL> pdf *= ( a - i ) * b / ( <NUM_LIT:1> - b ) / ( i + <NUM_LIT:1> ) <EOL> if pdf == <NUM_LIT:0.0> : break <EOL> cdf += pdf <EOL> cdf = min ( <NUM_LIT:1> , cdf ) <EOL> cdf = float ( "<STR_LIT>" % cdf ) <EOL> return cdf <EOL> else : <EOL> pdf = binomial_pdf ( x , a , b ) <EOL> cdf = pdf <EOL> for i in xrange ( x - <NUM_LIT:1> , - <NUM_LIT:1> , - <NUM_LIT:1> ) : <EOL> pdf /= ( a - i ) * b / ( <NUM_LIT:1> - b ) / ( i + <NUM_LIT:1> ) <EOL> if pdf == <NUM_LIT:0.0> : break <EOL> cdf += pdf <EOL> cdf = min ( <NUM_LIT:1> , cdf ) <EOL> cdf = float ( "<STR_LIT>" % cdf ) <EOL> return cdf <EOL> def binomial_cdf_inv ( cdf , a , b ) : <EOL> """<STR_LIT>""" <EOL> if cdf < <NUM_LIT:0> or cdf > <NUM_LIT:1> : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> cdf2 = <NUM_LIT:0> <EOL> a = int ( a ) <EOL> for x in xrange ( <NUM_LIT:0> , a + <NUM_LIT:1> ) : <EOL> pdf = binomial_pdf ( x , a , b ) <EOL> cdf2 = cdf2 + pdf <EOL> if cdf < cdf2 : <EOL> return x <EOL> return x <EOL> def binomial_pdf ( x , a , b ) : <EOL> """<STR_LIT>""" <EOL> if a < <NUM_LIT:1> : <EOL> return <NUM_LIT:0> <EOL> elif x < <NUM_LIT:0> or a < x : <EOL> return <NUM_LIT:0> <EOL> elif b == <NUM_LIT:0> : <EOL> if x == <NUM_LIT:0> : return <NUM_LIT:1> <EOL> else : return <NUM_LIT:0> <EOL> elif b == <NUM_LIT:1> : <EOL> if x == a : return <NUM_LIT:1> <EOL> else : return <NUM_LIT:0> <EOL> else : <EOL> if x > a - x : <EOL> p = <NUM_LIT:1> - b <EOL> mn = a - x <EOL> mx = x <EOL> else : <EOL> p = b <EOL> mn = x <EOL> mx = a - x <EOL> pdf = <NUM_LIT:1> <EOL> t = <NUM_LIT:0> <EOL> for q in xrange ( <NUM_LIT:1> , mn + <NUM_LIT:1> ) : <EOL> pdf *= ( a - q + <NUM_LIT:1> ) * p / ( mn - q + <NUM_LIT:1> ) <EOL> if pdf < <NUM_LIT> : <EOL> while pdf < <NUM_LIT> : <EOL> pdf /= <NUM_LIT:1> - p <EOL> t -= <NUM_LIT:1> <EOL> if pdf > <NUM_LIT> : <EOL> while pdf > <NUM_LIT> and t < mx : <EOL> pdf *= <NUM_LIT:1> - p <EOL> t += <NUM_LIT:1> <EOL> for i in xrange ( mx - t ) : <EOL> pdf *= <NUM_LIT:1> - p <EOL> pdf = float ( "<STR_LIT>" % pdf ) <EOL> return pdf <EOL> def facotrial ( n ) : <EOL> """<STR_LIT>""" <EOL> if n < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> fact = <NUM_LIT:1> <EOL> for i in xrange ( <NUM_LIT:2> , n + <NUM_LIT:1> ) : <EOL> fact = fact * i <EOL> return fact <EOL> def test ( ) : <EOL> print ( <NUM_LIT:1> - binomial_cdf ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) ) <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> test ( ) </s>
