Korrekte Anpassung mit scipy curve_fit einschließlich Fehler in x?

Question

Ich versuche ein Histogramm mit einigen Daten darin unter Verwendung scipy.optimize.curve_fit anzupassen. Wenn ich einen Fehler in y hinzufügen möchte, kann ich einfach einen weight auf die Anpassung anwenden. Aber wie kann man den Fehler in x anwenden (d. H. Der Fehler wegen Binning im Falle von Histogrammen)?

Meine Frage gilt auch für Fehler in x bei einer linearen Regression mit curve_fit oder polyfit; Ich weiß, wie man Fehler in y hinzufügt, aber nicht in x.

Hier ein Beispiel (zum Teil aus dem matplotlib documentation):

import numpy as np 
import pylab as P 
from scipy.optimize import curve_fit 

# create the data histogram 
mu, sigma = 200, 25 
x = mu + sigma*P.randn(10000) 

# define fit function 
def gauss(x, *p): 
    A, mu, sigma = p 
    return A*np.exp(-(x-mu)**2/(2*sigma**2)) 

# the histogram of the data 
n, bins, patches = P.hist(x, 50, histtype='step') 
sigma_n = np.sqrt(n) # Adding Poisson errors in y 
bin_centres = (bins[:-1] + bins[1:])/2 
sigma_x = (bins[1] - bins[0])/np.sqrt(12) # Binning error in x 
P.setp(patches, 'facecolor', 'g', 'alpha', 0.75) 

# fitting and plotting 
p0 = [700, 200, 25] 
popt, pcov = curve_fit(gauss, bin_centres, n, p0=p0, sigma=sigma_n, absolute_sigma=True) 
x = np.arange(100, 300, 0.5) 
fit = gauss(x, *popt) 
P.plot(x, fit, 'r--') 

Nun, diese passen (wenn er nicht ausfällt) hält die y-Fehler sigma_n, aber ich habe keine Möglichkeit gefunden mach es in Betracht ziehen sigma_x. Ich scannte ein paar Threads auf der scipy Mailingliste und fand heraus, wie man den absolute_sigma Wert und einen Beitrag auf Stackoverflow über asymmetrical errors verwendet, aber nichts über Fehler in beiden Richtungen. Ist es möglich zu erreichen?

Answer 1

scipy.optmize.curve_fit verwendet die standardmäßige nichtlineare Methode der kleinsten Quadrate und minimiert daher nur die Abweichung der Antwortvariablen. Wenn Sie einen Fehler in der zu berücksichtigenden unabhängigen Variable haben möchten, können Sie versuchen scipy.odr, die orthogonale Distanzregression verwendet. Wie sein Name andeutet, minimiert es sowohl unabhängige als auch abhängige Variablen.

Schauen Sie sich das folgende Beispiel an. Der fit_type Parameter bestimmt, ob scipy.odr volle ODR (fit_type=0) oder kleinste Quadrate Optimierung (fit_type=2) tut.

EDIT

Obwohl das Beispiel funktionierte es nicht viel Sinn hat, da die y-Daten auf den lauten x Daten berechnet wurde, die in einem Abstand ungleich indepenent Variable führte gerade. Ich aktualisierte das Beispiel, das jetzt auch zeigt, wie man RealData benutzt, das den Standardfehler der Daten anstelle der Gewichtungen vorgibt.

from scipy.odr import ODR, Model, Data, RealData 
import numpy as np 
from pylab import * 

def func(beta, x): 
    y = beta[0]+beta[1]*x+beta[2]*x**3 
    return y 

#generate data 
x = np.linspace(-3,2,100) 
y = func([-2.3,7.0,-4.0], x) 

# add some noise 
x += np.random.normal(scale=0.3, size=100) 
y += np.random.normal(scale=0.1, size=100) 

data = RealData(x, y, 0.3, 0.1) 
model = Model(func) 

odr = ODR(data, model, [1,0,0]) 
odr.set_job(fit_type=2) 
output = odr.run() 

xn = np.linspace(-3,2,50) 
yn = func(output.beta, xn) 
hold(True) 
plot(x,y,'ro') 
plot(xn,yn,'k-',label='leastsq') 
odr.set_job(fit_type=0) 
output = odr.run() 
yn = func(output.beta, xn) 
plot(xn,yn,'g-',label='odr') 
legend(loc=0)