Scipy und linearer Fit unter Beruecksichtigung von y-errors

[PaulchenPanther]

Hallo liebe Forumsgemeinde,

ich kaempfe nun schon seit zwei Tagen mit den Fitmoeglichkeiten in Python und hoffe nun hier auf Hilfe.
Mein Problem: Ich habe Messwerte, die grob eine Gerade darstellen. In y-Richtung gibt es einen Fehler, den ich beim Anfitten beruecksichtigen moechte.
Ich wuerde also gerne einen linearen Fit unter Beruecksichtigung des Fehlers meiner y-Werte vornehmen.

Ich lese schon eine ganze Weile in
1. http://wiki.scipy.org/Cookbook/LinearRegression
und
2. http://wiki.scipy.org/Cookbook/FittingData
Ich schaffe es aber nicht, die dort gefundenen Beispiele auf mein Problem zu uebertragen.
Beim ersten habe ich das Gefuehl, dass es eigentlich das ist, was ich moechte, aber leider wird "yerr" nicht beruecksichtigt.
Kann mir jemand einen Tipp in die richtige Richtung geben. Dankeschoen!

[EyDu]

Hallo,

was heißt den in deinem Fall "Fehler"? Gibt es da irgend eine Verteilung? Ansonsten werfe ich einfach mal Maximum Liklihood Estimation und Bayesian Linear Regression in den Raum.

[PaulchenPanther]

Die einzelnen y-Werte sind Means aus mehreren Einzelmessungen. Als dazugehoerigen Fehler nehme ich einfach die Standardabweichung.

[PaulchenPanther]

Ich glaube, ich habe gerade im Netz doch endlich ein Beispiel gefunden, bei dem genau das gemacht wird, was ich auch moechte.

Code: Alles auswählen

from pylab import *
from scipy.optimize import leastsq
from random import *
#Now we have to make a linear model. This function must take two arguments, the first one is a list of fit parameters and the second is a list of x values.
def func(p,x):
	return p[0]+p[1]*x
#Next we need to define the weighted residual of our fit. That is what is the difference between the data and model divided by the error bar on the data. 
#If you do not have error bars on your data this will just be the data minus the model.
def residuals(p,x,y,err):
	model=func(p,x)
	return (y-model)/err
#Now we will make a simulated data set. First we will make the true data, and then and some random noise to it. To show how much of a difference 
#the error bars make we will make every other data point have and error bar of 2 and the others with error bars of 10.
x=arange(10)
y=arange(10)
err=[2,10]*5
for i in range(len(err)):
	y[i]=y[i]+choice([-1,1])*random()*err[i]
#To use leastsq you first need to make a guess at the paramiters for your fit, we will call this p0. Then you call leastsq wiht the risidual function, p0, 
#args(other values that go into the residual funciton). This will then return a list of the best fit paramiters.
p0=[0,.5,0]
plsq=leastsq(residuals,p0,args=(x,y,err))
#Lets see how well we did. Remember the true values should be [0,1,0].
print plsq[0]
#As you can see the fit using the error bars did a much better job of finding the true values. We can also visualize this by plotting the data and fits on the same plot.
errorbar(x,y,yerr=err,fmt='o')
plot(x,func(plsq[0],x),'r-',label='Best fit within error bars')
plot(x,x,'b-',label='True fit')
legend(loc=2)
plt.show()

Dieses Beispiel werde dann mal auf meine eigenen Daten uebertragen.