Python, statistics: Test of proportions.

Our Python code implements the following:

description	Python code	Details
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
author: ernesto p. adorio
             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.