Pearson product-moment correlation

correlation.py

from math import sqrt
from random import sample
from datetime import datetime

def correlation(ser1, ser2):
    """Computes Pearson product-moment correlation coefficient for the given 2
    time series. See: http://bit.ly/1PHG1jy
    >>> from correlation import correlation
    >>> s1 = [random()*100 for i in range(100)]
    >>> _s1 = [100-s1[i] for i in range(100)]
    >>> assert abs(correlation(s1, s1) - 1.0) < 0.001
    >>> assert abs(correlation(s1, _s1) + 1.0) < 0.001
    """
    sum1 = sum2 = sum12 = sumsq1 = sumsq2 = 0
    n = len(ser1)
    for i, x1 in enumerate(ser1):
        global ali
        ali += 1
        x2 = ser2[i]
        sum1 += x1; sum2 += x2
        sum12 += x1*x2
        sumsq1 += x1**2; sumsq2 += x2**2
    nom = (n*sum12) - (sum1*sum2)
    den = sqrt(n*sumsq1 - sum1**2)*sqrt(n*sumsq2 - sum2**2)
    return nom/den

# Constants
NUM_SER = 200
SIZ_WIN = 9
STEPS = 12

def main():
    print 'Generating random data.'
    data = [sample(range(100), SIZ_WIN) for i in range(NUM_SER)]
    print 'Correlate.'
    time = datetime.now()
    for t in range(STEPS):
        R = []
        for i in range(NUM_SER):
            for j in range(i, NUM_SER):
                r = correlation(data[i], data[j])
                R.append((i, j, r))
        R.sort(key=lambda x: x[2])
        #print R[:10]
        for i, _ in enumerate(data):
            data[i] = sample(range(100), 1) + data[i][:-1]
    print 'Total correlation time: %s' % str(datetime.now()-time)

if __name__ == '__main__':
    main()