Untitled

#coding: utf8

'''
Compute Binomial distribution.

http://en.wikipedia.org/wiki/Binomial_distribution

times:

python2.7 28 secondes
pypy1.6nightly 1 minute 11 secondes

$ pypy1.6nightly --version
Python 2.7.1 (82bf0efcfe7d, Aug 11 2011, 02:06:34)
[PyPy 1.6.0-dev1 with GCC 4.4.3]
'''


from __future__ import division, print_function
import sys
from decimal import Decimal

class Factorial(object):
    '''
    Instanciate it, and use at as a factorial generator, but with cache.

    >>> fact = Factorial()
    >>> fact(4)
    24
    >>> fact(3) # should not compute anything
    6
    '''
    def __init__(self):
        self.factorials = [1, 1]

    def __call__(self, n):
        try:
            return self.factorials[n]
        except IndexError:
            current_index = len(self.factorials)

            while current_index <= n:
                self.factorials.append(current_index * self.factorials[-1])
                current_index += 1

            return self(n)

factorial = Factorial()

def C(n, k):
    '''
    Returns the combinason of *n*  *k*, C(n,k) = n! / ( k! (n-k)!).

    >>> # With a dice of three faces, how much combinaisons can I do with twos
    >>> # 12-13-23 == 3
    >>> C(3, 2)
    Decimal('3')
    '''
    # Use Decimal to avoid flooting point overflow
    return Decimal(factorial(int(n))) / Decimal(factorial(int(k)) * factorial(int(n-k)))

def P(n, k, p):
    """
    Repeat *n* time a random experiment with probability *p* of succes.

    P(n, k, p) returns the probability of getting *exactly* *k* sucessful
    experiment.

    >>> # what is the probability of getting only 1 six from a dice after 10 tries
    >>> 0 < P(10, 1, 1/6) < 0.4
    True
    """
    n, k, p = Decimal(n), Decimal(k), Decimal(p)

    return C(n, k) * p ** k * (1 - p) ** (n - k)

def P_sum(n, k, p):
    """
    Repeat *n* time a random experiment with probability *p* of succes.

    P_sum(n, k, p) returns the probability of getting *at least* *k* sucessful
    experiments.

    >>> # what is the probability of getting a six (or more) from a dice after 10 tries
    >>> 0.80 < P_sum(10, 1, 1/6) < 0.85
    True
    """
    return sum(P(n, k_p, p) for k_p in range(k, n + 1))

def optimise_max(k, p, min):
    '''
    Try to find the minimal number of try to get *at least* *k* sucessful event
    with a probability > min.

    The event itself has a probability of *p*.

    >>> # What is the minimal number of try I must do to be sure that I'll get \
        # a twos 6 (or more) with more than 50% chance.
    >>> optimise_max(2, 1/6, 0.5)
    10
    '''
    n = 1
    while True:
        P = P_sum(n, k, p)
        print(n, P, file=sys.stderr)
        if P > min:
            break
        n *= 2

    bound_a, bound_b = n // 2, n

    while bound_b - bound_a > 1:
        n = (bound_b + bound_a) // 2
        P = P_sum(n, k, p)

        print(n, P, file=sys.stderr)

        if P <= min:
            bound_a = n
        else:
            bound_b = n

    return n + 1

if __name__ == '__main__':
    print(optimise_max(1, 0.0009, 0.7))