Python, statistics: Ranking a single array of raw scores

There are some nonparametric statistical routines which expects rank values. If we are given raw scores
like X= [1,3,3,5,1,7], the rank values would be [1,2,2,3,1,4] if tied ranks get the same rank value.

But the ranks can be transformed to [0.5, 2.5, 2.5,3, .5, 4] if ties are replaced by the mean of tied ranks.

Here is Python code to compute ranks based on raw scores according to four strategies published in the reference.

"""
File      scores2ranks.py
Author    Ernesto Adorio, PhD.
               UPDEPP at Clarkfield, Pampanga
               ernesto.adorio@gmail.com
Desc      Conversion of raw scores to ranks using various strategies.
Version   0.0.1 October 1, 2011
License   Educational use only with proper attribution for research purposes.
Reference http://en.wikipedia.org/wiki/Ranking
"""

def scores2ranks(X, ztol = 1.0e-1, breakties = 1):
   """
    Converts raw scores to ranks, returning an array of ranks.
    Args
      X - scores to convert to ranks
      ztol - equality comparison tolerance
      breakties- strategy:
        0 - None.                               1234            ordinal ranking
        1 - replace  ties by mean of tied ranks.1 2.5, 2.5, 4 fractional ranking.
        2 - (competition rank)                  1224            standard competition ranking.
        3 - replace ties by highest tied rank.  1334            modified competition ranking. 
        4 - replace rank after ties in sequence 1223            dense ranking.
    References: For conversion of matrix scores to ranks using fractional ranking:
     http://my-other-life-as-programmer.blogspot.com/2011/02/python-converting-raw-scores-to-ranks.html         
   """
   Z = [(x, i) for i, x in enumerate(X)]
   Z.sort()
   n = len(Z)
   Rx = [0] * n 
   for j, (x,i) in enumerate(Z):
       Rx[i] = j+1
   if breakties == 0:
      return Rx
   s = 1           # sum of ties.
   start = end = 0 # starting and ending marks.
   for i in range(1, n):
       if abs(Z[i][0] -Z[i-1][0]) < ztol and i != n-1:
          pos = Z[i][1]
          s+= Rx[pos]
          end = i 
       else: #end of similar x values.
          if breakties == 1:
             tiedRank = float(s)/(end-start+1)
             for j in range(start, end+1):
                Rx[Z[j][1]] = tiedRank
          if breakties == 2 or  breakties == 4:
             tiedRank = Rx[Z[start][1]]      
             for j in range(start, end+1):
                Rx[Z[j][1]] = tiedRank
          if breakties == 3:
             tiedRank = Rx[Z[end][1]]      
          for j in range(start, end+1):
              Rx[Z[j][1]] = tiedRank
          start = end = i
          s = Rx[Z[i][1]]  
   
 
   if breakties == 4:
         #ensure that  the ranks are in sequence!
         for i, x in enumerate(sorted(list(set(Rx[:])))):
             for j, y in enumerate(Rx):
                 if y == x:
                    Rx[j] = i+1  
   return Rx
 

if __name__ == "__main__":
    X= [1,3,3,5,  1,  7]
    print "X = ", X
    print scores2ranks(X,  breakties = 4)

When the above code is run, it outputs

$ python scores2ranks.py 
X =  [1, 3, 3, 5, 1, 7]
strategy 0 : [1, 3, 4, 5, 2, 6]
strategy 1 : [1.5, 3.5, 3.5, 5.0, 1.5, 6]
strategy 2 : [1, 3, 3, 5, 1, 6]
strategy 3 : [2, 4, 4, 5, 2, 6]
strategy 4 : [1, 2, 2, 3, 1, 4]

I will be grateful if readers will discover any mistake.

Finite sampling demonstator

This is a mirror of an article in our extreme-solvers.blogspot.com

Given a finite population S with N elements which may not be unique (some elements are repeated),
we extract a finite sample X with n elements where n < N. The way we extract the n elements may be done in the following manner:

1. Permutions without replacement.
The ordering of the sample elements is important and a sample element once chosen in the sample may NOT be chosen again.
2. Permutations with replacement.
The ordering of the sample elements is important and a sample element once chosen in the sample may be chosen again.
3. Combinations without replacement
The ordering of the sample elements is NOT important and a sample element once chosen to be in the sample is not available again.
4. Combinations with replacement
The ordering of the sample elements is NOT important and an element chosen to be in the sample may still be chosen again.

The total number of samples for each type of finite sampling above is given in the following table:

Number of Samples

sampling	Combinations	Permutations
without replacement	$$\frac{N!}{n! (N-n)!}$$	$$\frac{N!}{(N-n)!}$$
with replacement	\frac{(N+n-1)}{n!(N-1)!}	$$N^n$$

Let $$S = [s_0, s_1, s_2, ....,s_{N-1}]$$. Our generated sample X is actually X = $$[s_{i_0}, s_{i_1}, ...., s_{i_{n-1}}]$$
where the indices $$i_0 to i_{n-1}$$ is sequentially generated by a combinatorial algorithm.

We may be interested in one the following statistic which is a random variable for the totality of all samples:

1. sample mean
2. sample sum or total
3. sample s.d.(standard deviation) (divisor is n-1)
4. populaton s.d(divisor is n)
5. sample var (sample variance)(divisor is n-1)
6. population variance(population variance) (divisor is n)
7. sample vaiance(sample variance) (divisor is n-1)
8. sample max (maximum value)
9. sample min (minimum value)
10. range (max - min )

We wish to gather data on ALL possible finite samples for the desire statistic,
the mean and standard deviation and of course the distribution table for the statistic which contain the
Columns for statistic, frequency, (rf) relative frequency, (crf) cumulative relative frequency. x rf

Here is a complete example for the sum of the numbers which show up in a throw of three dice:
The population consists of [1,2,3,4,5,6].
The population size is 6.
The sample size is 3.
The ordering is considered important, for example [1,3,2] will be considered different from [3,1,2]. Thus it
is a permutation with replacement.
The "first" sample is [1,1,1] with a total of 3 and the "last" sample is [6,6,6] with a total of 18.
To help with our computations, we use our solvers hosted at www.extreme. to do it for us!

We will only show the generated frequency distribution table as the 216 generated samples is too long for this page.

Sampling Statistic Frequency Distribution Table

x	f	rf	crf	x rf	(x-mu)^2 rf
3.0	1	0.00462962962963	0.00462962962963	0.0138888888889	0.260416666667
4.0	3	0.0138888888889	0.0185185185185	0.0555555555556	0.586805555556
5.0	6	0.0277777777778	0.0462962962963	0.138888888889	0.840277777778
6.0	10	0.0462962962963	0.0925925925926	0.277777777778	0.9375
7.0	15	0.0694444444444	0.162037037037	0.486111111111	0.850694444444
8.0	21	0.0972222222222	0.259259259259	0.777777777778	0.607638888889
9.0	25	0.115740740741	0.375	1.04166666667	0.260416666667
10.0	27	0.125	0.5	1.25	0.03125
11.0	27	0.125	0.625	1.375	0.03125
12.0	25	0.115740740741	0.740740740741	1.38888888889	0.260416666667
13.0	21	0.0972222222222	0.837962962963	1.26388888889	0.607638888889
14.0	15	0.0694444444444	0.907407407407	0.972222222222	0.850694444444
15.0	10	0.0462962962963	0.953703703704	0.694444444444	0.9375
16.0	6	0.0277777777778	0.981481481481	0.444444444444	0.840277777778
17.0	3	0.0138888888889	0.99537037037	0.236111111111	0.586805555556
18.0	1	0.00462962962963	1.0	0.0833333333333	0.260416666667
Sum	216	1.0		mean=10.5	variance=8.75
					std.dev=2.95803989155
Finite Population Parameters, Correction factor=0.774596669241
(N,n)	Mean	Pvar	Pstd	Svar	Sstd
(6, 3)	3.5	2.91666666667	1.70782512766	3.5	1.87082869339

Visit the solver Extreme Solvers: Stats sampling

Be sure to specify the right parameters for the solver for the above example, see the screen below:

Simple probability on cross-tab data problems.

I thought all of my ECON studes will be able to solve these problems but NO. Most of them expect canned packaged problems but are lead astray when they have to think more.

Consider the following frequency table of second year Engg students taking various subjects:

Subject\Student	Male	Female
Calculus	5	25
History	5	30
Physics	15	6

1. What is the probability that a 2nd Year Engg. student chosen is taking Calculus or is a Male?

2. What is the probability that a 2nd Year Enng. student is a Female or is taking History?

3. What is the probability that 2nd Year Engg. student is a female taking History?

4. What is the probability that a 2nd Year Engg. student is not taking Caculus?

Answers.
First form the full table with column and row sums.

Subject\Student	Male	Female	Sum
Calculus	15	25	40
History	5	30	35
Physics	15	6	21
Total	35	61	96

1.Use the fact that $$P(A\cup B) = P(A) + P(B) - P(A\capB)$$

The answer is then $$P(Calulus \cup Male) = P(Calculus) + P(Male) - P(Calculus \cap Male)$$
or $$\frac{40 + 35 - 15}{96}= \frac{60}{96}= 0.625$$

2. Similar to #1,$$P(History \cup Female) = P(Calculus) + P(Female) - P(History \cap Female)$$
Solving, $$\frac{61+35-30}{96}=\frac{66}{96}=\frac{11}{96}= 0.114583$$.

3. This is the simplest problem. From the above table, it is $$\frac{30}{96} = 0.3125$$

4. The complement is $$1 - P(Calculus)$$. Therefore $$1- \frac{40}{96} = \frac{54}{96} = 0.5625$$

The Poisson distribution as a limiting form of the Binomial Distribution

This entry is a mirror of the article in Digital Explorations. As you can see, the Wordpress Latex plugin works better compared to the Blogger plugin.


Simeon Dennis Poisson derived his distribution as a limiting form of the binomial distribution
as $$n \to \infty$$ while $$\lambda = n p $$ stays constant. It is very instructive for aspiring
mathematical statisticians to see the details.

We have $$p = \lambda/ n$$ and $$q = 1 - \lambda/n$$. The binomial distribution pmf
$$f(x)=\binom{n}{x} p^x q^{n-x}$$
then becomes $$\frac{n!}{x!(n-x)!} (\frac{\lambda}{n})^x (1 - \frac{\lambda}{n})^{n-x}$$ which when simplifying becomes $$\underbrace{\frac{n!}{x!(n-x)!}}_A \underbrace{(\frac{\lambda}{n})^x}_B
\underbrace{(1 - \frac{\lambda}{n})^{n-x}}_C$$.

Recall that the limit of a product is the
product of the limits of its factors. Consider C:
$$ C = (1 - \frac{\lambda}{n})^{n}(1 - \frac{\lambda}{n})^{-x}$$
As $$n\to \infty, (1 - \frac{\lambda}{n})^{n}$$ converges to $$e^{-\lambda}$$.
while $$(1 - \frac{\lambda/n})^{-x}$$ converges to 1.

Thus we have at the moment,
$$ A B e^{-\lambda}$$
Now consider $$A = \frac{n!}{x!(n-x)!}$$. This can be rewritten as
$$\frac{n(n-1) \cdots (n-k+1)(n-k)(n-k-1)\cdots 2 \cdot 1}{x! (n-x)!}$$ or
$$A= \frac{n(n-1) \cdots (n-(x+1))}{x!}$$ the numerator has x factors.
and factoring,
$$ A = \frac{n^k (1 (1 - 1/n) (1 - 2/n)\cdots (1 - (n-x+1)/n}{x!}
{x!}$$
The limit as $$n\to\infty$$ is thus $$A = \frac{n^x(1)1)\cdots(1)}{x!}$$
Putting together ABC, we now have

$$[\frac{n^x}{x!}][\frac{\lambda^n} {n^x}][e^{-\lambda}]$$

From which we have the pmf of the Poisson distribution
$$f(x) = \frac{ \lambda^x e^{-\lambda}}{x!}$$

Python, Statistics: The nonparametric Mann-Whitney Test

Draft! Untested. Do not use yet.

def mannwhitney(S1, S2):
    """
    Returns the Mann-Whitney U statistic of two samples S1 and S2.
    """
    # Form a single array with a categorical variable indicate the sample
    X = [(s, 0) for s in S1]
    X.extend([(s,1) for s in S2])
    R = Rank(X)

    # Compute needed parameters.
    n1 = len(S1)
    n2 = len(S2)

    # Compute total ranks for sample 1.          
    R1 = sum([R[i] for i, (x,j) in enumerate(X) if j == 0])
    u1 = R1 - (n1 + (n1+1)/2.0)
    u2 = n1 * n2 - u1
    U = min(u1, u2)

    mU     = n1 * n2 / 2.0
    sigmaU = sqrt((n1 *n2)*(n1 + n2 + 1)/12.0)
    return U, mu, sigmaU

Still needs to find resources for computing the discrete distribution function of the Mann-Whitney test. Blogger will appreciate any help. Failing this, the scipy module has a mannwhitneyu function which returns the U statistic and the p-value of the test.

Ranking an array

In reviewing nonparametric statistical procedures, we see that we are missing the Mann-Whitney test. Now scipy has a routine mannwhitneyu(x,y) which returns the U statistic and the pvalue of the test. The null hypothesis is that the two independent samples have the same medians.

It requires that the data is at least on the ordinal, interval or ratio scales or the data must be sortable or ranked. Since the values (after transformation) are ranks, we looked for our old routine for ranking and decided to revamp it. See Version 0.0.1

def rankarray(X, rankstartvalue = 1,  averageTies=True):
    """
    version 0.0.2 may 16, 2010
    """
    R = [rankstartvalue + i for i in range(len(X))]
    xi =  [(x,i) for i, x in enumerate(X)]
    xi.sort()
    if averageTies:
        start = 0
        end   = 0
        for i in range(1, len(X)):
           if xi[i][0] == xi[start][0]:
              end = i
           else:
              count = end-start + 1
              avgRank = xi[start][0] + xi[end][0]  / 2
              for j in range(start,end+1):
                  R[j] = avgRank
              start = i
              end   = i
              
        #Adjust for any trailing similar ranks.
        if start != end:
            count = end - start + 1
            for j in range(start,end+1):
                  R[j] = (R[start][0]+ R[end][0])/2
                  
    RR=[0] * len(X)
    for j,  (x,  i) in enumerate(xi):
        # print j, x, i,  R[i]
        RR[i] =R[j]
    return RR

if __name__ == "__main__":
    X = [7, 7, 4,4,4,4,4, 8,  6, 5, 1, 1]
    
    print "X=",  X
    print "Rank start value=%d, averageTies=%s"  %(0,True)
    R = rankarray(X,  rankstartvalue = 0,  averageTies = True)
    print R
    
    print "Rank start value=%d, averageTies=%s"  %(0,False)
    R = rankarray(X,  rankstartvalue = 0,  averageTies = False)
    print R
    
    print "Rank start value=%d, averageTies=%s"  %(1,True)
    R = rankarray(X,  rankstartvalue = 1,  averageTies = True)
    print R
    
    print "Rank start value=%d, averageTies=%s"  %(1,False)
    R = rankarray(X,  rankstartvalue = 1,  averageTies = False)
    
    print R
    

When the above script is run, it outputs:

Python 2.6.5 (r265:79063, Apr 16 2010, 13:57:41) 
[GCC 4.4.3] on toto-laptop, Standard
>>> X= [7, 7, 4, 4, 4, 4, 4, 8, 6, 5, 1, 1]
Rank start value=0, averageTies=True
[10, 10, 6, 6, 6, 6, 6, 11, 9, 7, 1, 1]
Rank start value=0, averageTies=False
[9, 10, 2, 3, 4, 5, 6, 11, 8, 7, 0, 1]
Rank start value=1, averageTies=True
[10, 10, 6, 6, 6, 6, 6, 12, 9, 7, 1, 1]
Rank start value=1, averageTies=False
[10, 11, 3, 4, 5, 6, 7, 12, 9, 8, 1, 2]

But we cannot guarantee 100 percent that there are no errors.

The recursive Durbin formula for the partial autocorrelation function.

This post is only for historical purposes and is superceded by a new direct CORRECTED implementation in Corrected post for partial autocorrelation function.This post may be deleted one week from now. Please visit the revised and corrected version.



Durbin developed a recursive formula for the pacfś given the acf($r_k$'s): \[ \phi_k= \phi_{k,k} = \begin{cases} r_1 & \text{for $k=1$}\\ \frac{r_k - \sum_{j=1}^{k-1} \phi_{k-1} r_{k-j} }{1 -\sum_{j=1}^{k-1}\phi_{k-1,j} r_j} & \text{for $k>1$}\\ \end{cases} \] \[ \[\phi_{k,j} = \begin{cases} \phi_{k-1, j} -\phi_{k}\phi_{k-1,k-j} & \text{for $j = 1,2, \ldots, k-1$}\\ \phi_{k,k} & \text{for $j =k$}\\ \end{cases} \]



Here is our initial Python implementation of the recursive function.

def durbin(r, lag = 1):
    #print"debug: inside durbin, lag=", lag
    
    def phi(k, j):
        #print "calling phi with ", k, j
        if k == j:
           if k == 1:
              return r[1]
           else:
              N= r[k] - sum([phi(k-1, k-1) * r[k-j] for j in range(1,k)])
              D= 1 - sum([phi(k-1, j) * r[j] for j in range(1, k)])
              return N/D
        else:
           return sum([phi(k-1,j) - phi(k,k) * phi(k-1,k-j) for j in range(1, k)])
    return phi(lag,lag) 

def  acf(X, k=1):
    """
    Computes the sample autocorrelation function coeffficient.
    for given lag k and input data X.
    """
    if k == 0:
        return 1.0
    flen = float(len(X))
    xbar = float(sum([x for x in X])) / flen
    D = sum([(x-xbar)**2 for x in X])
    N = sum([ (x-xbar)* (xtpk -xbar) for (x, xtpk) in zip(X[:-k],X[k:])])
    return N/D

def sacf(X,  maxlag):
    return [acf(X, i) for i in range(maxlag+1)] #includes 0 lag

#data from Danao's text. page 323.
X = ['772.9', '909.4', '1080.3', '1276.2', '1380.6', '1354.8', 
   '1096.9', '1066.7', '1108.7', '1109', '1203.7', '1328.2', '1380', 
   '1435.3', '1416.2', '1494.9', '1525.6', '1551.1', '1539.2', 
   '1629.1', '1665.3', '1708.7', '1799.4', '1873.3', '1973.3', 
   '2087.6', '2208.3', '2271.4', '2365.6', '2423.3', '2416.2', 
   '2484.8', '2608.5', '2744.1', '2729.3', '2695', '2826.7', 
   '2958.6', '3115.2', '3192.4', '3187.1', '3248.8', '3166', 
   '3279.1', '3489.9', '3585.2', '3676.5']

def Test1(maxlag):
    """
    Prints out values of acf and pacf for the Danao data set.
    """

    Y = [float(x) for x in X]
    r  = sacf(Y,  maxlag)
    for lag in range(1, 5):
      print r[lag], durbin(r, lag)
 
if __name__== "__main__":
   Test1(17)
 

The recursion is actually performed by an internal function phi() which accepts two arguments k, j. However, the output of the above even for low lag value of 5 is not so good( we may have miscoded the algorithm!). Anyway we present it in the hope that a gentle reader will show us the way to properly implement Durbin's algorithm! We will recheck and recheck the formulas.

A corrected version has been successfully written. Here is the wrong results returned by our first implementation.

$ python durbin.py
python durbin.py
0 0.925682317386 0.925682317386
1 0.852706579655 -0.0292160394675
2 0.787096604484 5.86791773408
3 0.737850083142 -3.27133094039

Notice the very bad value of 5.867+ for i = 2, lag = 3. No partial autocorrelations are ever greater than 1!

duh...

Python, Statistics: The non-parametric Runs test for randomness.

Top

1. Introduction
   
2. Probability mass function
   
3. Normal Approximation
   
4. Python Code
   
5. Swed and Eisenhart table
   

Top

Intro

Pmf

Python

Normal

Swed

Bottom

Introduction

Consider the sequence ABBAAABBABABBB. A subsequence of consecutive same letters is called a run. For our example, we use vertical bars to separate each run: A|BB|AAA|BB|A|B|A|BBB consists of 8 runs. This sequence contains 6 A's and 8 B's which we denote by (m,n)= (6,8). We would be suspicious of people claiming that the sequences AAAAAAABBBBBBB or ABABABABABABAB are random. The runs test was developed precisely to test the null hypothesis whether a sample of a sequence of binary symbols is randomly generated. The minimum number of runs of two symbols $m$ of one kind and $n$ of another kind is 2. On the other hand, the maximum possible number of runs is $2\, min(m, n)$ plus 1 if both m and n are different.

The probability mass function of the distribution of runs

The distribution of the number of runs has been "derived" or explained in [Mood, Graybill and Boes,"Introduction to Statistics" McGraw Hill, 3e, 1974] p.519-521. The probability mass function is given for any even number of runs $(z = 2k)$ as \[ prob(Z = z)=\frac{2\binom{m-1}{k-1}\binom{n-1}{k-1}}{\binom{m+n}{m}}\] and for an odd number of runs, $(z = 2k+1)$ as \[ prob(Z = z)=\frac{\binom{m-1}{k}\binom{n-1}{k-1} + \binom{m-1}{k-1} \binom{n-1}{k}}{\binom{m+n}{m}}\] Recall that $\binom{n}{r}$ is the number of combinations of $n$ things taken $r$ at a time.It can be computed as $\binom{n}{r} = \frac{n!}{r!(n-r)!}$ but this is inefficient, and instead we use the formula $\binom{n}{r} = \frac{n(n-1)\cdots (n-r+1)}{1\cdot 2\cdots r}$ It would be foolish to do by hand even with aid of a non-programmable scientific calculators, especialy when calculating the cumulative mass distribution function and is the reason we are writing Python codes.

Top

Intro

Pmf

Normal

Python

Swed

Bottom

The normal approximation to the distribution of runs

When m and n are both large, say m and n are both greater than 10, then an asymptotic normal approximation might be good enough. The mean is given by \[\mu = \frac{2mn}{m+n} +1\] and the variance would be given by \[\sigma^2 = \frac{2mn(2mn -m -n)}{(m+n)^2 (m+n-1}\]

Top

Intro

Pmf

Normal

Python

Swed

Bot

Python codes

Here is our latest version as of April 6. It is more complete now but as usual, users beware. It has not been thoroughly debugged.

Function	Description
nCr(n, r)	binomial coefficient
maxruns(m,n)	maximum runs is (m,n) sequence
pmfruns(z, m, n)	probability mass function
cmftableruns(m, n)	cumulative mass function table
pvalue(z, m, n, side=-1)	$p-$value of test statistic z
normrunstest(z, m, n, alpha, side)	normal approximation test
exactrunstest(z,m,n,alpha,side=0)	exact test for runs
swedeisenhart	Swed and Eisenhart (m,n)cumulative distribution table
isvalidruns(S)	return (m,n) if S is a valid runs sequence, else None

# -*- coding: utf-8 -*-
"""
file    testofruns.py
author  Ernesto P. Adorio
version 0.0.1 April 5, 2010
         
"""

from math import sqrt
from basictests import ztest,  ttest


def  nCr(n, r):
     """
     Computes the binomial coefficient, number of combinations of n things taken r at a time.
     """
     if r == 1:  return n
     if r == 0 or r == n: return 1
     
     if n < r   or n  < 0 or r  < 0: return 0
     if r > n-r: # find value of less looping.
        r = n-r
     f = 1
     for i in range(1, r+1):
        f *= (n - i+1) 
        f /= i
     return f
     
def maxruns(m,n):
     # maximum runs with m  symbols of the first kind and n symbols of the second.kind.
     return 2 * min(m, n) + (1 if m!=n else 0)
  
def pmfruns(z, m, n):
     """
     Computes the probability mass function for the distribution of runs.
     """
     if z < 2: 
        raise ValueError,  "in pmfruns, z must satisfy z >= 2"
     if  z % 2 ==0: #even z.
        k = z // 2 
        return 2.0 * nCr(m-1, k-1) * nCr (n-1, k-1)/ nCr(m+n, m)
     else : #odd z
        k = (z-1) // 2
        return float(nCr(m-1, k) * nCr(n-1, k-1) + nCr(m-1, k-1) * nCr(n-1, k)) / float(nCr(m + n, m))

def cmftableruns(m, n):
    """
    Returns a cumulative masss function for the distribution of runs.
    You must note that that the random variable z starts at 2 !
    """
    tot = 0
    result=[]
    
    maxz = maxruns(m,  n)
    for i in range(2, maxz+1):
       tot += pmfruns(i, m, n)
       result.append(tot)
    return  result


def pvalue(z, m,  n,  side=-1):
    "unchecked. to be tested first!"
    upperz = maxruns(m,n)
    if not (2 <= z <= maxruns) :
       raise ValueError, "in pvalue(), z must lie in (2, %d)" % upperz
    table = cmftableruns(m,n)
    if side == 0:
       p = table[z-2]
       if p > 0.5: 
          p = 1-0.5
       return p * 2
    elif side == -1:
        return table[z-2]
    elif side == 1:
       return 1 - table[z-2] 
        
def normrunstest(z,  m,  n,  alpha,  side):
    #Test using normal approximation, when both m and n are greater than 10.
    mean = 2.0 * m *n / (m+n) + 1
    variance = 2. * m * n * (2 *m *n -m -n)/((m+n)**2 * (m+n - 1))
    return ztest(z - mean,  sqrt(variance/n), alpha,  side) 


def swedeisenhart():
    print
    print """"""      for i in range(2, 11):         print '' % i,     print "
"     for m in range(2, 11):         for n in range(m, 11):             rowvals = cmftableruns(m, n)             maxcol = len(rowvals)             print "
" % (m,n),             for j in range(min(maxcol, 9)):                 print "" % rowvals[j],             for j in range(maxcol, 9):                 print "",             print "
"     print "
P(Z <= a) when H_0 is true in the Runs test

(m, n)
%d
(%d,%d)
%5.3f
"

def exactrunstest(z, m, n, alpha, side=0):
    #returns pvalue, test stat z and lower and upper limits.
    cmf = cmftableruns(m,n)
    print cmf, 
    a = b = None
    if side ==0:
       alpha /= 2.0
       alpha1 = alpha
       alpha2 = 1-alpha1
    elif side == 1:
       alpha2 = 1-alpha
    elif side == 0:
       alpha1 = alpha
    for i in range(len(cmf)):
       if cmf[i] <= alpha1:
          a  = i+2 
       if cmf[i] >= alpha2:
          b =  i+2
    p = pvalue(z,  m,n, side)
    if side == 0:
       return (p, z, (a,b))
    elif side == 1:
       return (p, z, b)
    elif side == -1:
       return (p, z, a)

   
def isvalidruns(S):
        D = {}
        for s in S:
           if s in D:
                D[s] +=1
           else:
                D[s] = 1
                if len(D) != 2:
                   raise  ValueError, "in isvalidruns(): more than two symbols"
        
        counts = [] 
        for (i,s) in enumerate(D):
             counts.append(D[s])
        a,b  = counts[0], counts[1]
        return min(a, b), max(a,b)
 
  
  
 
def Test():
   swedeisenhart() # decomment line to turn on table generation test.

   #print "Exact runs test for z = 5, (m, n)= (5,7)"
   #print exactrunstest(5, 5,7, 0.05, side=0)
   """
   print pmfruns(3, 5,4)
   print normrunstest(5, 4,9,0.05, -1)
   print "Table for runs test."
   m, n = 3,5
   print cmftableruns(m,n)
   """

if __name__ == "__main__":
   Test()

The Python module above calls the ztest which was already described in previous posts. We will add all files required in one zip file for convenience later.

Top

Intro

Pmf

Python

Normal

Swed

Bot

Swed and Eisenhart Table

Here is an html table obtained by running swedeisenhart() routine.

$P(Z <= a)$ when $H_0$ is true in the Runs test

(m, n)	2	3	4	5	6	7	8	9	10
(2,2)	0.333	0.667	1.000
(2,3)	0.200	0.500	0.900	1.000
(2,4)	0.133	0.400	0.800	1.000
(2,5)	0.095	0.333	0.714	1.000
(2,6)	0.071	0.286	0.643	1.000
(2,7)	0.056	0.250	0.583	1.000
(2,8)	0.044	0.222	0.533	1.000
(2,9)	0.036	0.200	0.491	1.000
(2,10)	0.030	0.182	0.455	1.000
(3,3)	0.100	0.300	0.700	0.900	1.000
(3,4)	0.057	0.200	0.543	0.800	0.971	1.000
(3,5)	0.036	0.143	0.429	0.714	0.929	1.000
(3,6)	0.024	0.107	0.345	0.643	0.881	1.000
(3,7)	0.017	0.083	0.283	0.583	0.833	1.000
(3,8)	0.012	0.067	0.236	0.533	0.788	1.000
(3,9)	0.009	0.055	0.200	0.491	0.745	1.000
(3,10)	0.007	0.045	0.171	0.455	0.706	1.000
(4,4)	0.029	0.114	0.371	0.629	0.886	0.971	1.000
(4,5)	0.016	0.071	0.262	0.500	0.786	0.929	0.992	1.000
(4,6)	0.010	0.048	0.190	0.405	0.690	0.881	0.976	1.000
(4,7)	0.006	0.033	0.142	0.333	0.606	0.833	0.955	1.000
(4,8)	0.004	0.024	0.109	0.279	0.533	0.788	0.929	1.000
(4,9)	0.003	0.018	0.085	0.236	0.471	0.745	0.902	1.000
(4,10)	0.002	0.014	0.068	0.203	0.419	0.706	0.874	1.000
(5,5)	0.008	0.040	0.167	0.357	0.643	0.833	0.960	0.992	1.000
(5,6)	0.004	0.024	0.110	0.262	0.522	0.738	0.911	0.976	0.998
(5,7)	0.003	0.015	0.076	0.197	0.424	0.652	0.854	0.955	0.992
(5,8)	0.002	0.010	0.054	0.152	0.347	0.576	0.793	0.929	0.984
(5,9)	0.001	0.007	0.039	0.119	0.287	0.510	0.734	0.902	0.972
(5,10)	0.001	0.005	0.029	0.095	0.239	0.455	0.678	0.874	0.958
(6,6)	0.002	0.013	0.067	0.175	0.392	0.608	0.825	0.933	0.987
(6,7)	0.001	0.008	0.043	0.121	0.296	0.500	0.733	0.879	0.966
(6,8)	0.001	0.005	0.028	0.086	0.226	0.413	0.646	0.821	0.937
(6,9)	0.000	0.003	0.019	0.063	0.175	0.343	0.566	0.762	0.902
(6,10)	0.000	0.002	0.013	0.047	0.137	0.287	0.497	0.706	0.864
(7,7)	0.001	0.004	0.025	0.078	0.209	0.383	0.617	0.791	0.922
(7,8)	0.000	0.002	0.015	0.051	0.149	0.296	0.514	0.704	0.867
(7,9)	0.000	0.001	0.010	0.035	0.108	0.231	0.427	0.622	0.806
(7,10)	0.000	0.001	0.006	0.024	0.080	0.182	0.355	0.549	0.743
(8,8)	0.000	0.001	0.009	0.032	0.100	0.214	0.405	0.595	0.786
(8,9)	0.000	0.001	0.005	0.020	0.069	0.157	0.319	0.500	0.702
(8,10)	0.000	0.000	0.003	0.013	0.048	0.117	0.251	0.419	0.621
(9,9)	0.000	0.000	0.003	0.012	0.044	0.109	0.238	0.399	0.601
(9,10)	0.000	0.000	0.002	0.008	0.029	0.077	0.179	0.319	0.510
(10,10)	0.000	0.000	0.001	0.004	0.019	0.051	0.128	0.242	0.414

If you have any suggestions, or corrections, or features request do not hesitate to email the author at ernesto.adorio @ gmail . com.

Bottom

Top

Intro

Pmf

Python

Normal

Swed

Bot

I have encountered the html problem: Your HTML cannot be accepted: Closing tag has no matching opening tag: A ???!!! I wonder where its !!

Testing the significance of the correlation coefficient r

The sample correlation coefficient is a measure of the strength of the linear relationship between two variables X and Y. Actually it is the coefficient of determination $r^2$ which is a better measure. It is computed by the formula

$$ r =\frac{\sum_{i = 1}^n (x_i - \overline{x} ) (y_i - \overline{y})}{S_x S_y} $$

where $S_x$ and $S_y$ are the standard deviations of the X and Y samples. We have For $S_x = \sqrt{\frac{\sum (x_i -\overline{x})^2}{n-1}}$ with a similar formula for $S_y$. The value of $r$ varies from -1 (perfect negative correlation) to +1 (perfect positive correlation). A negative correlation means that as one variable increases, the other other variable decreases, or in other words, both variables varies in opposite directions. On the other hand a positive correlation arises when both variables varies together in the same direction: if one variable increases, the other also increases and if one variable decreases, the other also decreases.

Correlation is generally used when it is difficult to determine the causative order. Does X causes Y or does Y causes X? It is heavily used in the Psychology field for example. In this post, we describe how to test the significance of the computed correlation coefficient.

The sample test statistic $t= \frac{r-\rho}{\sqrt{\frac{1-r^2}{n-2}}}$ follow assymptotically the t-distribution with $df=n-2$ degree of freedom and tables of the t-distribution may be used. Thus to write our Python program to use the ttest routine described previously in our blog, we specify the sample test statistic and the standard error.

def corrtest( rho, r, n, alpha = 0.5, side = 0):
     stattest = (r - rho)
     se = sqrt( (1-r*r)/(n-2.0))
     return ttest(r-rho, se, alpha, side)

Neat, is not it? We obtain the p-value, the distribution test statistic and the critical limits(values). For side = 0, a double sided test, the two limits do not strictly form the true confidence interval. Here is the complete testforcorr.py Python file where we included the ttest function. It calls the rstats.py module we have written before to use R-names for the Python statistical functions.

# -*- coding: utf-8 -*-
"""
file     testofcorr.py
author   Ernesto P. Adorio
         ernesto.adorio @ gmail.com
         UP at Clarkfield
         Angeles, Pampanga

version  0.0.1 april 4, 2010
"""

from math import *
from rstats import *



def ttest(samplestat, se, df, alpha= 0.05, side=0):
    """                                            
    T-test of a sample statistic.                  
    Arguments:                                     
     samplestat- sample statistic                  
     se -    standard error of sample statistic    
     df      degree of freedom                     
     alpha - significance level                    
     side  - -1 left-sided test                    
             0  double sided test                  
             +1 right-sided test                   
    """                                            
    Ttest = samplestat/se                          
    if side ==0:                                   
        pvalue = pt(Ttest, df)                     
        if Ttest > 0:                              
            pvalue = 1-pvalue                      
        pvalue *=2                                 
        tcrit1 =    qt(alpha/2,   df)              
        tcrit2 =    qt(1-alpha/2.0,df)             
        return pvalue, Ttest, (tcrit1,tcrit2)      
    elif side == -1:                               
       pvalue = pt(Ttest,  df)                     
       tcrit = qt(alpha,  df)                      
       return pvalue, Ttest, tcrit                 
    else:                                          
       pvalue = 1- pt(Ttest,  df)                  
       tcrit  = qt(1.0-alpha,  df)                 
       return pvalue, Ttest, tcrit                 


def cov(X,Y):
    n = len(X)
    if n != len(Y):
        raise "ArgumentError", "in cov: len(X) != len(Y)"
    xbar = float(sum(X))/ n
    ybar = float(sum(Y))/n
    return sum([(x -xbar)*(y-ybar) for x,y in zip(X,Y)])/(n-1)

def Sx(X):
    xbar = float(sum([x for x in X]))/n
    return sqrt( sum([(x - xbar)**2 for x in X]) /(n-1.0))
 
def cor(X,Y): 
    # correlation coefficient of X and Y. 
    return cov(X,Y)/(Sx(X)*Sx(Y))


def corrtest( rho, r, n, alpha = 0.5, side = 0):
     stattest = (r - rho)
     se = sqrt( (1-r*r)/(n-2.0))
     return ttest(r-rho, se, n-2,  alpha, side)


def xycorrtest(X,Y, rho, alpha = 0.5, side= 0):
     if len(X) != len(Y):
        raise ArgumentError, "in xycorrtest: X and Y must have same length"
     r = cor(X,Y)
     return corrtest(rho, r, len(X), alpha, side)


def Test():
     #Berenson , p. 546 
     rho = 0
     r = 0.951
     alpha = 0.05
     side = 1
     n = 14
     print "Berenson's example, p. 546"
     print corrtest( rho, r, n, alpha, side)

if __name__ == "__main__":
    Test()

When we invoke the Python interpreter it will run the Berenson, p. 546 example. The output is


toto@toto-laptop:~/Blogs/statistics$ python testofcorr.py
Berenson's example, p. 546
(8.9865653363219167e-08, 10.654779470653486, 1.7822875556491591)


The extremely low p-value means that we can reject the Null hypothesis that the correlation coefficient is zero.
To support this decision, the obtained distribution test statistic is greater then 1.78 critical value limit, i.e.,
is in the rejection region.


I need to specify a darker background for <code> blocks and a central download repository for
the files. Stay tuned. We will add more useful technical contents soon!

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.

t-test for means

Whenever the variance of the normal population is unknown, then the t-test is appropriate for testing the Null hypothesis that the population mean is some value $\mu_0$. The variance itself is estimated from the finite sample:

$$ s^2 = \frac{\sum_{i=1}^n (x_i - \overline{x})^2}{n-1}$$

The null hypothesis : $$ H_0 = \mu = \mu_0$$
The alternative hypostheses (select only 1!) are :$$\mu < \mu, \mu \ne \mu_0, \mu > \mu_0$$
The test statistic : $$t_{test}= \frac{\overline{x} - \mu_0}{\sqrt{s^2 /n}}$$

The critical values are found by consulting the t distribution tables with an additional parameter df for degrees of freedom or simply computed.

The following Python function ttestcomp performs the necessary computations given the hypothesized population mean, the significance level alpha, the kind of test.

It should return the pvalue, the test statistic and the critical value(s). For a two sided kind of test,the critical values form a pair, otherwise only one critical value is returned.

from math import *

from math import *
from rstats import *


def ttestcomp(X, popmean0, alpha= 0.05, side= 0):
    """
    Arguments:
     X  - input array
     popmean0- assumed population mean, 
     alpha - significance level
     side  -1 - left sided test
            0 - double sided test
            1 - right sided test
    Return value:
     p-value, ztest, critical values.
    """
    n = len(X)
    xbar = sum(X)/ float(n)
    s2   = sum([(x - xbar)**2  for x in X])/float(n)
     
    ttest = (xbar - popmean0) / sqrt(s2/(n-1.0))
    df = n -1      

    if side ==0:
       pvalue = pt(ttest, df)/2.0
       tcrit1 = qt(alpha/2, df) 
       tcrit2 = qt(1-alpha/2.0, df)
       return pvalue, ttest, (tcrit1,tcrit2)
    elif side == -1:
       pvalue = pt(ttest, df)
       tcrit  = qt(alpha)
       return pvalue, ttest, tcrit
    else:
       pvalue = 1- pt(ttest, df)
       tcrit  = qt(1-alpha,df) 
       return pvalue, ttest, tcrit

if __name__ == "__main__":
   popmean0 = 7
   popvar  = 36
   samplesize = 20
   alpha = 0.05
   side = 0
   X= rnorm(10, sqrt(popvar), size = samplesize) # note the difference in arguments with R rnomrm

   print "t-test for means"
   print "finished generating X"
   print X
   print "two-sided:", ttestcomp(X, popmean0, alpha, side)
   print "right-sided:", ttestcomp(X, popmean0, alpha, side=1)

The test code included in our Python code performs the ttest for both two-sided and right-sided alternatives. It first generates a random sample of size 20 from a normal distribution with mean 7 and variance 36.

When the above program is run, the output is

python ttestcomp.py
t-test for means
finished generating X
[ 4.7861463 17.75238035 9.82834098 20.48616187 16.22260671
20.04664846 6.02676326 18.18174814 10.20050244 -3.87442282
2.58404396 10.68370551 10.27102987 15.81674596 17.14011521
18.65564402 18.28997982 9.6119956 7.5671266 11.56864725]
two-sided: (0.49936270026320678, 3.4725570605535814, (-2.0930240544082634, 2.093024054408263))
right-sided: (0.0012745994735864352, 3.4725570605535814, 1.7291328115213671)

Note that for the two-sided test, we accept the Null hypothesis that the population mean is 7 whereas for the right-sided test, we REJECT the Null hypothesis that the population mean is 7. This is not unexpected since we generated a random sample with mean 10.

The rstats.py may be obtained from our Wordpress blog at rstats.py.


We hasten to add that the t-test may already be available in scipy.stats but we present our own for pedagogical purposes.

Z Test for means.

I was looking at the index to statistics when I discovered there were no codes for basic hypothesis testing!

Now the mean of a sample of size $n$ has best estimate value or test statistic $\overline{x} = \sum x_i$ with a standard error of $\srt{\sigma / n}$ where $\sigma^2$ is the known population variance. If the variance of the population is unknown, then it can be estimated from the values of the sample.

\[ s^2 = \frac{\sum_{i=1}^n (x_i - \overline{x})^2}{n-1} \]

The assumption is that the sample is obtained from a normal population with variance known.

The null hypothesis : $H_0 = \mu = \mu_0$
The alternative hypostheses (select only 1!) are :$\mu < \mu, \mu \ne \mu_0, \mu > \mu_0$
The test statistic : $Z_{test} \frac{\overline{x} - \mu_0}{\sqrt{\sigma^2 /n}}$

The critical values are found by consulting the normal distribution tables or simply computed.
For a two-sided test, the test alpha ($\alpha$) is half that of the specified value.


The following Python code performs the necessary computations given the hypothesized population mean,
the population variance, the significance level alpha, the kind of test and the test statistic.
The test function performs the ztest for both two-sided and right-sided alternatives.

It should return the pvalue, the test statistic and the critical value(s). For a two sided kind of test,the critical values form a pair, otherwise only one critical value is returned.

from math import *

from rstats import *
def ztestcomp(popmean0, popvar, samplesize, xbar, alpha= 0.05, side= 0):
    """
    Arguments:
     xbar  - sample mean
     popmean0- assumed population mean, 
     popvar - population variance
     alpha - significance level.
     side  -1 - left sided test
            0 - double sided test
            1 - right sided test
     xbar - sample mean
    Return value:
     p-value, ztest, critical values.
    """
    ztest = (xbar - popmean0) / sqrt(popvar/(samplesize-1.0))
          
    if side ==0:
       pvalue = pnorm(abs(ztest))/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-alpha) 
       return pvalue, ztest, zcrit

if __name__ == "__main__":
   popmean0 = 5
   popvar  = 36
   samplesize = 20
   alpha = 0.10
   side = 0
   xbar = 7

   print "two-sided:", ztestcomp(popmean0, popvar, samplesize, xbar, alpha, side)
  
   print "right-sided:", ztestcomp(popmean0, popvar, samplesize, xbar, alpha, side=1)
   
  

When the above program is run, the output is

~/Blogs/statistics$ python tests.py

toto@toto-laptop:~/Blogs/statistics$ python tests.py
two-sided: (0.42067237303427146, 1.0, (-1.9599639845400545, 1.959963984540054))
right-sided: (0.15865525393145707, 1.0, 1.6448536269514722)

Since the $p$-value is greater than the specified significance level (5 %= 0.05), we accept, or "do not reject" the Null hypothesis that the true population mean is 5.

Textbooks recommend the ztest only for n > 30. Our sample size is smaller, and is only used for illustration.

The rstats.py may be obtained from our Wordpress blog at rstats.py.

The ztestcomp() above only performs the necessary computations. A separate routine should be written if an input input sample array is desired, and is left as an exercise.

We hasten to add that the z-test may already be available in scipy.stats but we present our own for pedagogical purposes.

Test your Statistics Understanding! Part I

This test is mirrored in the Digital Explorations weblog. The latex plugin for Blogger is still not up to par with those available for Wordpress.This test is mirrored in Dr-Adorio-Adventures I don't really know what's wrong with this Latex plugin.It prefers you build up the equation, render it, re-edit, rerender.... What works in Wordpress may not work at once in Blogger!

From the Word Bank provided at the end of the problems, write down the best answer.

1. The value computed by $\frac{(\sum x-\overline{x})(\sum y - \overline{y})}{(n-1)\sigma_x \sigma_y}$
2. A plot which quickly shows the distribution of the data, showing minimum, Q1, median, Q2, maximum values.
3. A plot which simply shows points $(x_i, y_i)$ to help discern patterns between a variable $X$ and $Y$.
4. In a sample, the value with the greatest frequency.
5. The value $(x-\mu)/\sigma$ corresponding to a raw value $x$.
6. In ordinary regression, the value computed by $\overline{y}- b \overline{x}$.
7. The probability of commiting a type I error.
8. The probability of obtaining at least an extreme value as the test statistic observed.
9. A statement that there is no difference between parameters of two populations.
10. In ordinary regression, the value computed by $\frac{n\sum xy -\sum x\sum y}{n\sum x^2 -(\sum x)^2}$
11. The type of error commited when failing to reject the null hypothesis when in fact it is False.
12. A french mathematician who lived in England who was the first to apply the normal curve equation.
13. A bell shaped, symmetric continuous distribution also called gaussian distribution.
14. The discrete distirbution which describes the probability of success of a single event.
15. A continuous bell shaped distribution which arises when the standard deviation of the population is not known.
16. The discrete distribution which describes rare events like deaths due to kick of horses.
17. The distribution which describes the number of having babies until a boy appears.

Most of the answers are found in the following Word Bank:

p-value, significance level, scatterplot, boxplot, correlation coefficient, 
mode,  z-score, intercept, beta coefficient, slope, Null hypothesis, 
Type II error, degree of freedom,Legendre, de Moivre,  
Normal distribution curve, t-distribution, Bernoulli distribution, 
Poisson distribution, geometric distribuition

Answers to be published on April 8, 2010 in this blog entry page.

Care to comment? The latex renderer looks ok now...