# Untitled

By: a guest on Aug 11th, 2011  |  syntax: Python  |  size: 3.12 KB  |  hits: 63  |  expires: Never
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
1. #coding: utf8
2.
3. '''
4. Compute Binomial distribution.
5.
6. http://en.wikipedia.org/wiki/Binomial_distribution
7.
8. times:
9.
10. python2.7 28 secondes
11. pypy1.6nightly 1 minute 11 secondes
12.
13. \$ pypy1.6nightly --version
14. Python 2.7.1 (82bf0efcfe7d, Aug 11 2011, 02:06:34)
15. [PyPy 1.6.0-dev1 with GCC 4.4.3]
16. '''
17.
18.
19. from __future__ import division, print_function
20. import sys
21. from decimal import Decimal
22.
23. class Factorial(object):
24.     '''
25.    Instanciate it, and use at as a factorial generator, but with cache.
26.
27.    >>> fact = Factorial()
28.    >>> fact(4)
29.    24
30.    >>> fact(3) # should not compute anything
31.    6
32.    '''
33.     def __init__(self):
34.         self.factorials = [1, 1]
35.
36.     def __call__(self, n):
37.         try:
38.             return self.factorials[n]
39.         except IndexError:
40.             current_index = len(self.factorials)
41.
42.             while current_index <= n:
43.                 self.factorials.append(current_index * self.factorials[-1])
44.                 current_index += 1
45.
46.             return self(n)
47.
48. factorial = Factorial()
49.
50. def C(n, k):
51.     '''
52.    Returns the combinason of *n*  *k*, C(n,k) = n! / ( k! (n-k)!).
53.
54.    >>> # With a dice of three faces, how much combinaisons can I do with twos
55.    >>> # 12-13-23 == 3
56.    >>> C(3, 2)
57.    Decimal('3')
58.    '''
59.     # Use Decimal to avoid flooting point overflow
60.     return Decimal(factorial(int(n))) / Decimal(factorial(int(k)) * factorial(int(n-k)))
61.
62. def P(n, k, p):
63.     """
64.    Repeat *n* time a random experiment with probability *p* of succes.
65.
66.    P(n, k, p) returns the probability of getting *exactly* *k* sucessful
67.    experiment.
68.
69.    >>> # what is the probability of getting only 1 six from a dice after 10 tries
70.    >>> 0 < P(10, 1, 1/6) < 0.4
71.    True
72.    """
73.     n, k, p = Decimal(n), Decimal(k), Decimal(p)
74.
75.     return C(n, k) * p ** k * (1 - p) ** (n - k)
76.
77. def P_sum(n, k, p):
78.     """
79.    Repeat *n* time a random experiment with probability *p* of succes.
80.
81.    P_sum(n, k, p) returns the probability of getting *at least* *k* sucessful
82.    experiments.
83.
84.    >>> # what is the probability of getting a six (or more) from a dice after 10 tries
85.    >>> 0.80 < P_sum(10, 1, 1/6) < 0.85
86.    True
87.    """
88.     return sum(P(n, k_p, p) for k_p in range(k, n + 1))
89.
90. def optimise_max(k, p, min):
91.     '''
92.    Try to find the minimal number of try to get *at least* *k* sucessful event
93.    with a probability > min.
94.
95.    The event itself has a probability of *p*.
96.
97.    >>> # What is the minimal number of try I must do to be sure that I'll get \
98.        # a twos 6 (or more) with more than 50% chance.
99.    >>> optimise_max(2, 1/6, 0.5)
100.    10
101.    '''
102.     n = 1
103.     while True:
104.         P = P_sum(n, k, p)
105.         print(n, P, file=sys.stderr)
106.         if P > min:
107.             break
108.         n *= 2
109.
110.     bound_a, bound_b = n // 2, n
111.
112.     while bound_b - bound_a > 1:
113.         n = (bound_b + bound_a) // 2
114.         P = P_sum(n, k, p)
115.
116.         print(n, P, file=sys.stderr)
117.
118.         if P <= min:
119.             bound_a = n
120.         else:
121.             bound_b = n
122.
123.     return n + 1
124.
125. if __name__ == '__main__':
126.     print(optimise_max(1, 0.0009, 0.7))