Lineare Anpassung in Python mit Unsicherheit in x- und y-Koordinaten

Question

Hallo Ich möchte meine Kollegen Python-Benutzer fragen, wie sie ihre lineare Anpassung durchführen.

Ich habe seit den letzten zwei Wochen auf Methoden/Bibliotheken suchen, diese Aufgabe zu erfüllen, und ich möchte meine Erfahrung teilen:

Wenn Sie eine lineare Anpassung auf der Methode der kleinsten Quadrate basierend ausführen möchten Sie habe viele Möglichkeiten. Zum Beispiel können Sie Klassen in numpy und scipy finden. Mich I haben durch den einen präsentiert von linanp (die folgt die Gestaltung der linanp Funktion in IDL) entschieden:

http://nbviewer.ipython.org/github/djpine/linfit/blob/master/linfit.ipynb

Dieses Verfahren nimmt man die Sigmas im y-Achse sind die Einführung koordiniert Ihre Daten anzupassen .

Wenn Sie jedoch die Unsicherheit in der x- und y-Achse quantifiziert haben, gibt es nicht so viele Optionen. (Es gibt kein IDL "Fitexy" Äquivalent in den wissenschaftlichen Hauptbibliotheken von Python). Bisher habe ich nur die "kmpfit" -Bibliothek gefunden, um diese Aufgabe zu erfüllen. Glücklicherweise hat es eine sehr komplette Website beschreibt alle seine Funktionalität:

https://github.com/josephmeiring/kmpfit http://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#

Wenn jemand weitere Ansätze weiß, dass ich lieben würde, sich als gut zu kennen.

In jedem Fall hoffe ich, dass dies hilft.

Answer 1

Orthogonal distance regression in Scipy ermöglicht Ihnen die nichtlineare Anpassung mit Fehlern in x und y.

Unten sehen Sie ein einfaches Beispiel, das auf dem Beispiel auf der scipy-Seite basiert. Es versucht, einigen randomisierten Daten eine quadratische Funktion zuzuordnen.

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.odr import * 

import random 

# Initiate some data, giving some randomness using random.random(). 
x = np.array([0, 1, 2, 3, 4, 5]) 
y = np.array([i**2 + random.random() for i in x]) 

x_err = np.array([random.random() for i in x]) 
y_err = np.array([random.random() for i in x]) 

# Define a function (quadratic in our case) to fit the data with. 
def quad_func(p, x): 
    m, c = p 
    return m*x**2 + c 

# Create a model for fitting. 
quad_model = Model(quad_func) 

# Create a RealData object using our initiated data from above. 
data = RealData(x, y, sx=x_err, sy=y_err) 

# Set up ODR with the model and data. 
odr = ODR(data, quad_model, beta0=[0., 1.]) 

# Run the regression. 
out = odr.run() 

# Use the in-built pprint method to give us results. 
out.pprint() 
'''Beta: [ 1.01781493 0.48498006] 
Beta Std Error: [ 0.00390799 0.03660941] 
Beta Covariance: [[ 0.00241322 -0.01420883] 
[-0.01420883 0.21177597]] 
Residual Variance: 0.00632861634898189 
Inverse Condition #: 0.4195196193536024 
Reason(s) for Halting: 
    Sum of squares convergence''' 

x_fit = np.linspace(x[0], x[-1], 1000) 
y_fit = quad_func(out.beta, x_fit) 

plt.errorbar(x, y, xerr=x_err, yerr=y_err, linestyle='None', marker='x') 
plt.plot(x_fit, y_fit) 

plt.show() 

Answer 2

Sie können mit dem größten Eigenwert Eigenvektor der Kovarianzmatrix verwenden, um lineare Anpassung auszuführen.

import numpy as np 
import matplotlib.pyplot as plt 

x = np.arange(6, dtype=float) 
y = 3*x + 2 
x += np.random.randn(6)/10 
y += np.random.randn(6)/10 

xm = x.mean() 
ym = y.mean() 

C = np.cov([x-xm,y-ym]) 
evals,evecs = np.linalg.eig(C) 

a = evecs[1,evals.argmax()]/evecs[0,evals.argmax()] 
b = ym-a*xm 

xx=np.linspace(0,5,100) 
yy=a*xx+b 

plt.plot(x,y,'ro',xx,yy) 
plt.show()