-
Notifications
You must be signed in to change notification settings - Fork 0
/
binomial_test.py
83 lines (67 loc) · 2.43 KB
/
binomial_test.py
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
"""A simple implementation of the Binomial test.
Copyright Emanuele Olivetti, 2011.
This program is distributed under the GNU General Public Licence v3.0.
"""
import numpy as np
import scipy.misc as sm
import matplotlib.pyplot as plt
class Binomial(object):
"""Binomial distribution.
PMF, CDF and test accept both scalars and vectors as input.
This implementation is necessary since scipy.stats.binom as of
v0.7.0 does not work.
"""
def __init__(self, n, p):
self.n = n
self.p = p
def __call__(self, k):
"""Returns PMF value in k.
"""
return self.pmf(k)
def pmf(self, k):
"""Binomial PMF evaluate in k.
"""
return sm.comb(self.n, k) * self.p**k * (1.0 - self.p)**(self.n-k)
def cdf(self, k):
"""Binomial CDF in k.
"""
k = np.atleast_1d(np.round(k))
return np.array([self.pmf(np.arange(ki+1)).sum() for ki in k])
def test(self, k):
"""Binomial test, i.e., P(k<=e) under Binomial distribution.
"""
return self.cdf(k)
if __name__=='__main__':
testset_size = 50
epsilon = 0.5
num_errors = 18
alpha = 0.05
# not working in scipy 0.7.0:
# import scipy.stats as ss
# binomial = ss.binom(n=testset_size, pr=epsilon)
# binomial.pdf(num_errors)
# This motivates the implementation of Binomial in this file.
print "Binomial Test:"
binomial = Binomial(n=testset_size, p=epsilon)
p_value = binomial.test(num_errors)
print "P(k <= e=%s | testset_size=%s, epsilon=%s) = %s" % (num_errors, testset_size, epsilon, p_value)
print "Confidence level is alpha =", alpha
print "The Binomial test",
if p_value <= alpha :
print "REJECTS the null hypothesis of no information."
else:
print "DOES NOT REJECT the hypothesis of no information."
k = np.arange(testset_size+1)
plt.figure()
plt.plot(k, binomial.pmf(k), 'o')
plt.xlabel('number of errors (e)')
plt.ylabel('Binomial(%s, %s) PMF' % (testset_size, epsilon))
plt.figure()
plt.plot(k, binomial.cdf(k), 'o')
plt.xlabel('number of errors (e)')
plt.ylabel('Binomial(%s, %s) CDF' % (testset_size, epsilon))
plt.figure()
plt.plot(k,binomial.test(k), 'o')
plt.title("Binomial test, testset_size=%i epsilon=%s" % (testset_size, epsilon))
plt.xlabel('number of errors (e)')
plt.ylabel('p(k,<=e) under Binomial(%s, %s)' % (testset_size, epsilon))