ich finde es manchmal nervig, daß es in Python keine fertige Routine gibt, die eine lineare Regression so rechnet wie ich es will: Mit Berücksichtigung der Fehlerfortpflanzung und Angabe für Fehler von Steigung und Achsenabschnitt. (Jede Korrektur dieser Aussage in Form eines Links ist herzlich willkommen
.)Zumindest für den zweiten Teil des Problems habe ich mal eine schnelle Lösung zusammengehackt:
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()Viel Spaß damit,
Christian
edit: Kleine Korrektur (7.3.07) - Funktionalität voll erhalten und nicht geändert.