Das deutsche Python-Forum

Code: Alles auswählen

from numpy import arange
from math import sqrt
import pylab

#prepare data
x = arange(0.0, 2.0, 0.05)
nse = 0.3 * pylab.randn(len(x))
y = 2 + 3 * x + nse

#do the fit without considering errors using pylab.polyfit (same as in scipy)
m, b = pylab.polyfit(x, y, 1)

#plot data and polyfit fit
pylab.plot(x, y, 'bo', x, m * x + b, 'k-', linewidth = 2)

#caculate sums of squares:
n = len(x) #assuming len(x) == len(y)
xmean = x.mean()
ymean = y.mean()
ssxx = (x**2).sum() - n * xmean**2
ssyy = (y**2).sum() - n * ymean**2
ssxy = (x*y).sum() - n * xmean * ymean

#calculate intercept
m2 = ssxy / ssxx
#calculate slope (a)
b2 = ymean - m2 * xmean
#calculate linear regression coefficient
r2 = ssxy**2 / (ssxx * ssyy)

#calculate estimator for the variance:
s = sqrt(abs((ssyy - (ssxy**2 / ssxx)) / n - 2))
#calculate standard errors of intercept and slope
seb = s * sqrt((1/n) + (xmean**2 / ssxx))
sem = s / sqrt(ssxx)

#print report:
print "polyfits intercept: %.3f and slope: %.3f" % (m, b)
print "full intercept: %.3f +/- %.3f and slope: %.3f+/-%.3f and coefficient: %.3f" % (m2, sem, b2, seb, r2)

#plot all
pylab.plot(x, m2 * x + b2, 'r-', linewidth = 1)

pylab.show()