Check ordinary least squares calculation for correctness

Question

I've written ordinary least squares fitting of market data for learning purposes, so performance is not an issue. I'd like a reviewer to look at correctness and programming flaws, if any, in the code:

#!/usr/bin/env python

import math
import sys
import csv
import traceback

import ipdb

def excepthook_pm(typ, value, tb):
    traceback.print_exception(typ, value, tb)
    ipdb.pm()

sys.excepthook = excepthook_pm


class ExpandingStdDev(object):

    """
        Pure Python implementation of http://www.johndcook.com/standard_deviation.html

        It seems to be accurate enough, but it's totally inefficient since it's pure
        Python and uses dot access everywhere (costly in Python).
    """

    def __init__(self, xvals=[]):
        self.k = 0
        self.M_k_prev = 0.0
        self.M_k = 0.0
        self.S_k_prev = 0.0
        self.S_k = 0.0
        self.xvals = []
        if xvals:
            self.add_xvals(xvals)

    def push(self, x):
        self.xvals.append(x)
        self.k += 1
        if self.k == 1:
            self.M_k_prev = self.M_k = x
            self.S_k_prev = 0.0
        else:
            self.M_k = self.M_k_prev + (x - self.M_k_prev) / self.k
            self.S_k = self.S_k_prev + (x - self.M_k_prev) * (x - self.M_k)
            self.M_k_prev = self.M_k
            self.S_k_prev = self.S_k

    def add_xvals(self, xvals):
        for x in xvals:
            self.push(x)

    def mean(self):
        if self.k > 0:
            return self.M_k
        return 0.0

    def variance(self):
        if self.k > 1:
            return self.S_k / (self.k - 1)
        return 0.0

    def stddev(self):
        return math.sqrt(self.variance())


class OLS(object):

    def __init__(self, xvals, yvals):
        self.xvals = xvals
        self.yvals = yvals
        self.slope = None
        self.intercept = None

    def calculate(self):
        """ Basic ordinary least squares calculation. """
        xvals, yvals = self.xvals, self.yvals
        sumx, sumy = map(sum, [xvals, yvals])
        sumxy = sum(map(lambda x, y: x*y, xvals, yvals))
        sumxsq = sum(map(lambda x: x**2, xvals))
        Nsamp = len(xvals)
        # y = a*x + b
        # a (slope)
        slope = (Nsamp*sumxy - sumx*sumy) / (Nsamp*sumxsq - sumx**2)
        # b (intercept)
        intercept = (sumy - slope*sumx) / Nsamp
        self.slope = slope
        self.intercept = intercept
        return slope, intercept

    def predict(self, x):
        return self.slope*x + self.intercept


class TrendLinesBasic(object):

    def __init__(self, yvals, yvals_open, window, delimiter=';'):
        self.yvals = yvals
        self.yvals_open = yvals_open
        self.window = window
        self.delimiter = delimiter

    def find_all_trendpoints(self, fname):
        with open(fname, 'wb') as fo:
            w = csv.writer(fo, delimiter=self.delimiter)
            for idx in range(0, len(self.yvals), self.window):
                yvals_cut = self.yvals_open[idx:idx+self.window-1]
                ols = OLS(range(idx, idx+self.window-1), yvals_cut)
                ols.calculate()
                sd = ExpandingStdDev(yvals_cut).stddev()
                trend_start = ' '.join(self.yvals[idx][:2])
                trend_fin = ' '.join(self.yvals[idx+self.window-1][:2])
                y0 = ols.predict(idx)
                yn = ols.predict(idx + self.window-1)
                line = [trend_start, y0, y0+2*sd, y0-2*sd, trend_fin, yn, yn+2*sd, yn-2*sd]
                print 'Line:', ';'.join(map(str, line))
                w.writerow(line)

def read_data(fpath):
    with open(fpath, 'rb') as fo:
        rdr = csv.reader(fo)
        lines = [x for x in rdr]
        return lines


if __name__ == '__main__':
    # file columns: 'date', 'hour', 'open', 'high', 'low', 'close', 'volume'
    datfile = 'market_data.csv'
    yvals = read_data(datfile)
    yvals_open = map(float, [y[-2] for y in yvals])
    _N = len(yvals)
    _n = 10
    window = _N/_n
    print 'N: {0} n: {1} window: {2}'.format(_N, _n, window)
    tlines = TrendLinesBasic(yvals, yvals_open, window)
    outfname = '{0}_trendlines.csv'.format(datfile.replace('.csv', ''))
    tlines.find_all_trendpoints(outfname)
    print 'Wrote file:', outfname

Answer 1

Here are a few pointers:

1. Start function/method names with a verb. So I would change mean and variance to get_mean and get_variance respectively.

2. Your return statements in mean and variance can be put into one line:

   def mean(self):
       return 0.0 if self.k <= 0 else self.M_k
   

3. Your names could be more descriptive. This is an implementation of mathematical equations and unfortunately mathematic naming conventions don't translate well to code. Try to find better, more descriptive variable names.

4. Spacing is your friend. You have several places where you spacing is good. However, when you get into calculations, you bunch everything together. Try and separate logical sections by blank lines.

As my final point, simply look over PEP8, Python's official style guide. It'll help style your code better than I could.