Skip to content

Instantly share code, notes, and snippets.

@riccardoscalco
Last active July 22, 2022 16:32
Show Gist options
  • Save riccardoscalco/5356167 to your computer and use it in GitHub Desktop.
Save riccardoscalco/5356167 to your computer and use it in GitHub Desktop.
Confidence intervals on linear regression

Confidence intervals on linear regression

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

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
@Guitlle
Copy link

Guitlle commented Feb 12, 2016

Great. Thanks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment