Confidence intervals on linear regression
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 NonePython code for the evaluation of linear regression and confidence intervals between two random variables x and y.