Testování hypotéz - Binomický test

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