probably-probability/Binomial.py at main · DolevArtzi/probably-probability · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
from RandomVariable import RandomVariable
from Uniform import Uniform
import math

"""
Binomial RV

Number of successes in n independent trails with success probability p in each trial
"""
class Binomial(RandomVariable):
    def __init__(self, n, p,symmetric = True):
        super().__init__()
        self.min = 0
        self.max = n
        self.n = n
        self.p = p
        self.qn = math.pow(1 - p, n)
        self.params.append(n)
        self.params.append(p)
        self.name = 'binomial'
        if p == .5:
            self.setSymmetric(True)

    def mgf(self):
        return lambda t: (1-self.p + self.p * math.e ** t) ** self.n

    def pdf(self, k):
        if k < 0 or k > self.n or k % 1:
            return 0
        return math.comb(self.n, k) * math.pow(self.p, k) * math.pow(1 - self.p, self.n - k)

    def expectedValue(self):
        return self.n * self.p
    def _expectedValue(self,*params):
        print(params)
        print(self)
        n = params[0]
        p = params[1]
        return n*p

    def _valid(self,*params):
        n = params[0]
        p = params[1]
        return 0 <= p <= 1 and 0 <= n == int(n)

    def variance(self):
        return self.n * self.p * (1 - self.p)

    """
    Generates a Bin(n,p) random variable using the Inverse Transform Method
    (PnC book pg. 231, "Simulation" [Ross] pg. 57).

    Exploits the following recursive property for pdfs of binomials:

        p_x(0) = (1-p)^n

        p_x(i+1) = p/[1-p] * [n-i]/[i+1] * p_x(i)

        * Currently left as an artifact, not meant for use. Use genVar instead.
    """
    def _binomialInverseTransform(self):
        C = self.p / (1 - self.p)
        U = Uniform(0, 1).genVar()
        pr = self.qn
        F = pr
        i = 0
        while U >= F:
            pr *= (C * (self.n - i) / (i + 1))
            F += pr
            if i == self.n:
                break
            i += 1
        return i

    def _f(self,i):
        return (self.n - i) / (i+1)
    def genVar(self):
        return self.inverseTransform(self.p / (1 - self.p), self.qn , self._f)

    """Uses same recursive trick as above"""
    def cdf(self, k):
        if k < 0:
            return 0
        if k > self.n:
            return 1
        C = self.p / (1 - self.p)
        pr = self.qn
        F = pr
        i = 0
        while i < k:
            pr *= (C * (self.n - i) / (i + 1))
            F += pr
            i += 1
        return min(F, 1)

    def _getMomentRelatedFunction(self):
        return self.mgf()

if __name__ == '__main__':
    x = Binomial(20,0.5)
    print(x.cdf(15))