## Friday, April 2, 2010

### Python, statistics: Test of proportions.

Our Python code implements the following:

description Python codeDetails
cum. distribution function pbinom(x, n, p)$P(X \le x);f(x)= \binom{n}{x}p^xq^{n-x)}$
critical value onesidedcritbinom(n,p,alpha) $P(X\le xrit) \le \alpha$
critical values twosidedcritbinom(n,p,alpha)$P(xcrit1\le X\le xcrit2) \le \alpha$
exact test of proportion exactproptest(p0, x, n, alpha, side = 0) test using binomial pdf
approx test of prop.normproptest(p0, x, n, alpha, side = 0) test using normal approximation.
comparison of two prop. twosampleproptest(x1,n1, x2,n2,alpha,side)two sample test of prop.
test cases Walpole(),Walpole_normaltests() Test cases from "Intro. to statistics,3e" book.

"""
file:      testofproportions.py
UP Clarkfield
Pampanga
version 0.0.1 april, 3, 2010
"""
from math import  *
from rstats import  *
import scipy.stats as stat

def ztest(samplestat, se, alpha =0.05, side=0):
"""
Normal test of a sample statistic.
Arguments:
teststat- teststatistic
se -    standard error of sample statistic
alpha - significance level
side  - -1 left-sided test
0  double sided test
+1 right-sided test
"""
Ztest = samplestat/se
print "Ztest=",  Ztest
if side ==0:
pvalue = pnorm(Ztest)
if Ztest > 0.0:
pvalue = 1.0 -pvalue
pvalue *= 2.0
zcrit1 = qnorm(alpha/2)
zcrit2 = qnorm(1-alpha/2.0)
return pvalue, Ztest, (zcrit1,zcrit2)
elif side == -1:
pvalue = pnorm(Ztest)
zcrit = qnorm(alpha)
return pvalue, Ztest, zcrit
else:
pvalue = 1- pnorm(Ztest)
zcrit  = qnorm(1.0-alpha)
return pvalue, Ztest, zcrit

def pbinom(x,  n,  p):
"""
Computest the cumulative probability density function
of the binomial distribution up :  P(X <= x)
"""
#print "input to pbinom:",  x,  n,  p
q = 1.0 - p
pdf = cdf = q ** n
f = p/q
for i in range (1,  x+1):
pdf *= ((n -i + 1.0)/i * f)
cdf += pdf
return cdf

def onesidedcritbinom(n,  p,  alpha):
"""
Determines critical value xcrit such that P(X <= xcrit) < alpha
"""
q = 1.0 - p
pdf = cdf = q ** n
f = p/q
xcrit = None
for i in range (1,  x+1):
if xcrit is None and cdf >= alpha:
xcrit = i - 1
pdf *= ((n -i + 1.0)/i * f)
cdf += pdf
return xcrit

def twosidedcritbinom(n,  p,  alpha):
"""
Determines value xrit1, xrit2 such that P(X < xcrit1) = alpha/2 andP(X>xcrit2) = alpha/2
"""
if  not (0 <= alpha <= 1.0):
return None
q = 1.0 - p
pdf = cdf = q ** n
f = p/q
xcrit1 = None
xcrit2 = None
alpha1 = alpha/2.0
alpha2 = 1.0- alpha1
for i in range (1,  n+1):
if xcrit1 is None and cdf > alpha1:
print cdf,  alpha1
xcrit1 = i - 1
if xcrit2 is None and cdf > alpha2:
xcrit2 = i-1
break
pdf *= ((n -i + 1.0)/i * f)
cdf += pdf
return (xcrit1,  xcrit2)

def exactproptest(p0, x, n, alpha, side = 0):
"""
p0- assumed population proportion
x  - number of success of desired characteristic in sample
n  - sample size
Returns pvalue, teststatistic, critical value
"""
p = float(x) / n
if side == 0:
pvalue = pbinom(x,  n,  p0)
if pvalue > 0.5:
pvalue = 1 - pvalue
return pvalue,  x,  twosidedcritbinom(n, p0,   alpha)

elif side == -1:
pvalue = pbinom(x,  n,  p0)
return pvalue,   x,  onesidedcritbinom(n,  p0,  alpha)
elif side == 1:
pvalue = 1-0,  pbinom(x,  n,  1.0-p0)
return pvalue,  x,  onesidedcritbinom(n,  p0,  alpha)

def normproptest(p0,  x, n, alpha,  side = 0):
print "inputargs:",  p0,  x,  n,  alpha,  side
samplestat = x - n * p0
se = sqrt(n * p0 * (1-p0))
print "normproptest:", samplestat, "se=", samplestat,  se
return ztest(samplestat,  se,  alpha,  side)

def twosampleproptest(x1,  n1,   x2,  n2,  alpha,  side):
p1hat = float(x1)/n1
p2hat = float(x2)/n2
phat = float(x1 +x2)/(n1+n2)
print "p1hat, p2hat",  p1hat,  p2hat
samplestat =  p1hat - p2hat
se = sqrt(phat*(1.0-phat) * (1.0/n1 + 1.0/n2))
return ztest(samplestat,  se,  alpha,  side)

def Walpole_tests():
print "Example 10, p.326"
p0  = 0.7
alpha = 0.10
side = 0
x = 8
n = 15
return exactproptest(p0,  x,  n,  alpha,  side)

def Walpole_normaltests():
print "Example 11, p.329"
p0  = 0.6
alpha = 0.05
side = 1
x = 70
n = 100
print normproptest(p0,  x,  n,  alpha,  side)
print
print "Example12,  p. 330"
x1,  n1 = 120,  200
x2,  n2 = 240,  500
alpha = 0.025
side = 1
print twosampleproptest(x1,  n1,  x2,  n2,  alpha,  side)

if __name__ == "__main__":
print 'testing pbinom'
print Walpole_tests()
print Walpole_normaltests()


When Python(via the Eric Python IDE) runs the script file above, it outputs the following:

Python 2.6.4 (r264:75706, Dec  7 2009, 18:43:55)
[GCC 4.4.1] on toto-laptop, Standard
>>> testing pbinom
Example 10, p.326
0.0500125400538 0.05
(0.13114257338312119, 8, (7, 13))
Example 11, p.329
inputargs: 0.6 70 100 0.05 1
normproptest: 10.0 se= 10.0 4.89897948557
Ztest= 2.04124145232
(0.020613416668581852, 2.0412414523193152, 1.6448536269514722)

Example12,  p. 330
p1hat, p2hat 0.6 0.48
Ztest= 2.86972021592
(0.0020541757121497195, 2.8697202159177571, 1.959963984540054)
None


We will explain further in the days ahead how to test proportions as we are at the moment focusing on writing the codes.

#### 1 comment:

1. This comment has been removed by the author.