Máme čtyři události 0, 1, A, B. 0 a 1 se navzájem vylučují. Stejně tak A a B. Cíl je otestovat hypotézu, zda pravděpodobnost A za podmínky 1 je větší, než pravděpodobnost A za podmínky 0. Jinými slovy: Zda událost 1 zvyšuje pravděpodobnost události A. Matematicky v tom není problém, numericky ano, protože se tam vyskytují velká čísla. A někdy víc než jen velká. V Pythonu však není problém.
Program má na vstupu 5 čísel:
1) číslo "n0a" udávající počet výskytů události A za předpokladu události 0
2) číslo "n0b" udávající počet výskytů události B za předpokladu události 0
3) číslo "n1a" udávající počet výskytů události A za předpokladu události 1
4) číslo "n1b" udávající počet výskytů události B za předpokladu události 1
5) práh pravděpodobnosti pro zamítnutí HA
Standardní volba HA je 0.01 nebo 0.05. Volba 0.5 je daleko za hranicí, která odděluje odvahu a šílenství. Volba nad 0.5 je zcela nesmyslná, stejně jako pod 0.0.
# -*- coding: cp1250 -*-
import sys
import math
def div(a, b):
r = 0.0
e = 1.0
for i in range(100):
d = a / b
r = r + e * d
a = 10 * (a - d * b)
e = e / 10.0
return r
def C(n, k):
return math.factorial(n) / (math.factorial(k) * math.factorial(n - k))
def I01(ae):
sa = 0
sb = 1
for (a, e) in ae:
b = e + 1
sa = b * sa + a * sb
sb = b * sb
return (sa, sb)
def gen_ae(na, nb):
ae = []
for k in range(0, nb + 1):
ae.append((C(nb, k) * ((-1) ** k), na + k))
return ae
def suma(ab):
sa = 0
sb = 1
for (a, b) in ab:
sa = b * sa + a * sb
sb = b * sb
return (sa, sb)
n0a = int(sys.argv[1])
n0b = int(sys.argv[2])
n1a = int(sys.argv[3])
n1b = int(sys.argv[4])
alpha = float(sys.argv[5])
print
print "n0a =", n0a, "n0b =", n0b
print "n1a =", n1a, "n1b =", n1b
print "alpha =", alpha
print
ae0 = gen_ae(n0a, n0b)
ae1 = gen_ae(n1a, n1b)
(I0a, I0b) = I01(ae0)
(I1a, I1b) = I01(ae1)
ab = []
for (a0, e0) in ae0:
for (a1, e1) in ae1:
ab.append((a0 * a1, (e0 + 1) * (e0 + e1 + 2)))
(Ia, Ib) = suma(ab)
pA = 1 - div(Ia * I0b * I1b, Ib * I0a * I1a)
print
print "H0: p(A|1) > p(A|0)"
print "HA: p(A|1) <= p(A|0)"
print
print "p(HA) =", pA
print
if pA <= alpha:
print u"Doporučujeme zamítnou HA a přijmout H0."
if pA >= 1 - alpha:
print u"Doporučujeme zamítnou H0 a přijmout HA."
if pA > alpha and pA < 1 - alpha:
print u"nelze rozhodnout mezi H0 a HA."
print