Confidence intervals on linear regression

README.md

Confidence intervals on linear regression

Python code for the evaluation of linear regression and confidence intervals between two random variables x and y.

linreg.py

import scipy.stats, numpy
def linear_regression(x, y, prob):
    """
    Return the linear regression parameters and their <prob> confidence intervals.

    ex:
    >>> linear_regression([.1,.2,.3],[10,11,11.5],0.95)
    """
    x = numpy.array(x)
    y = numpy.array(y)
    n = len(x)
    xy = x * y
    xx = x * x

    # estimates

    b1 = (xy.mean() - x.mean() * y.mean()) / (xx.mean() - x.mean()**2)
    b0 = y.mean() - b1 * x.mean()
    s2 = 1./n * sum([(y[i] - b0 - b1 * x[i])**2 for i in xrange(n)])
    print 'b0 = ',b0
    print 'b1 = ',b1
    print 's2 = ',s2
    
    #confidence intervals
    
    alpha = 1 - prob
    c1 = scipy.stats.chi2.ppf(alpha/2.,n-2)
    c2 = scipy.stats.chi2.ppf(1-alpha/2.,n-2)
    print 'the confidence interval of s2 is: ',[n*s2/c2,n*s2/c1]
    
    c = -1 * scipy.stats.t.ppf(alpha/2.,n-2)
    bb1 = c * (s2 / ((n-2) * (xx.mean() - (x.mean())**2)))**.5
    print 'the confidence interval of b1 is: ',[b1-bb1,b1+bb1]
    
    bb0 = c * ((s2 / (n-2)) * (1 + (x.mean())**2 / (xx.mean() - (x.mean())**2)))**.5
    print 'the confidence interval of b0 is: ',[b0-bb0,b0+bb0]
    return None

Great. Thanks.