úterý 15. března 2011

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



2 komentáře:

al-Quaknaa řekl(a)...

Nerozumím úplně přesně tomu, jak ten program funguje (především proto, že jsem tomu nedal moc času), ale jako pythonista si říkám, že by ho značně zjednodušilo využití některých funkcí z math.factorial (http://docs.python.org/library/math.html), a**n místo pow(a, n) a a/b místo div(a, b). Ale to je tak stupidní, že pro to asi máte nějaké důvody, kterým akorát nerozumím...?

Taky by se hodil indenting, ale to je zase nejspíš chyba blogger.com.

Giovanni řekl(a)...

math.factorial a ** jsem nevyužil, protože jsem měl obavu, aby fungoval pro velká čísla jako třeba 100000. Zbytečně. Vyzkoušel jsem a funguje. Neceločíselné dělení pro velká čísla opravdu nefunguje. Namísto div by šlo použít float(fractions.Fraction(numerator, denominator)), ale nechám to už tak.