Lineare Anpassung einschließlich aller Fehler mit NumPy / SciPy

Eine Methode, die gut funktioniert und tatsächlich ein besseres Ergebnis liefert, ist die Bootstrap-Methode. Wenn Datenpunkte mit Fehlern angegeben werden, verwendet man einen parametrischen Bootstrap, und jeder x und y Wert beschreiben eine Gaußsche Verteilung. Dann zeichnet man einen Punkt aus jeder dieser Distributionen und erhält ein neues Bootstrapped-Sample. Die Ausführung einer einfachen nicht gewichteten Anpassung ergibt einen Wert für die Parameter.

Dieser Vorgang wird einige 300 bis einige tausend Mal wiederholt. Man wird am Ende eine Verteilung der Fit-Parameter erhalten, bei der man Mittelwert und Standardabweichung nehmen kann, um Wert und Fehler zu erhalten.

Eine weitere nette Sache ist, dass man dadurch keine einzige Fitkurve erhält, aber viele davon. Für jeden interpolierten x -Wert kann man wieder Mittelwert und Standardabweichung der vielen Werte f(x, param) nehmen und ein Fehlerband erhalten:

Weitere Schritte in der Analyse werden dann noch hunderte Male mit den verschiedenen Anpassungsparametern durchgeführt. Dies wird dann auch die Korrelation der Fit-Parameter berücksichtigen, wie man in der obigen Abbildung deutlich sehen kann: Obwohl eine symmetrische Funktion an die Daten angepasst wurde, ist das Fehlerband asymmetrisch. Dies bedeutet, dass interpolierte Werte auf der linken Seite eine größere Unsicherheit aufweisen als auf der rechten Seite.

Bitte beachten Sie, dass aus der Dokumentation von curvefit :

  

sigma: Keine oder N-Länge-Sequenz       Wenn nicht, wird dieser Vektor als relative Gewichtung in der       Kleinste-Quadrate-Problem.

Der entscheidende Punkt hier ist als relative Gewichtung , daher sollten yerr in Zeile 53 und 2*yerr in 57 Ihnen ähnliche, wenn nicht das gleiche Ergebnis geben.

Wenn Sie den tatsächlich Restfehler erhöhen, werden die Werte in der Kovarianzmatrix größer. Sagen wir, wenn wir die y += random in y += 5*random in der Funktion generate_data() :

Vergleicht das ursprüngliche Ergebnis:

Beachten Sie auch, dass die Parameterschätzung nun weiter von (2,3) abweicht, wie wir es von einem erhöhten Restfehler und einem größeren Konfidenzintervall von Parameterschätzungen erwarten würden.

Kurze Antwort

Für absolute Werte, die die Unsicherheit in y enthalten (und in x für den Fall odr):

• Verwenden Sie in der scipy.odr -Frage stddev = numpy.sqrt(numpy.diag(cov)) wo die cov ist die Kovarianzmatrix odr gibt in der Ausgabe.
• Verwenden Sie in der scipy.optimize.curve_fit -Frage absolute_sigma=True
  Flagge.

Für relative Werte (schließt Ungewissheit aus):

• Verwenden Sie im Fall scipy.odr den sd-Wert aus der Ausgabe.

• Verwenden Sie im Fall scipy.optimize.curve_fit absolute_sigma=False flag.

• Verwenden Sie numpy.polyfit wie folgt:

p, cov = numpy.polyfit(x, y, 1,cov = True) errorbars = numpy.sqrt(numpy.diag(cov))

Lange Antwort

Es gibt etwas undokumentiertes Verhalten in allen Funktionen. Meine Vermutung ist, dass die Funktionen relative und absolute Werte mischen. Am Ende ist diese Antwort der Code, der entweder angibt, was Sie wollen (oder nicht), basierend darauf, wie Sie die Ausgabe verarbeiten (es gibt einen Fehler?). Außerdem könnte curve_fit kürzlich das Flag 'absolute_sigma' erhalten haben?

Mein Punkt ist in der Ausgabe. Es scheint, dass odr die Standardabweichung berechnet, da es keine Unsicherheiten gibt, ähnlich wie bei Polyfit, aber wenn die Standardabweichung aus der Kovarianzmatrix berechnet wird, sind die Unsicherheiten vorhanden. Die curve_fit tut dies mit absolute_sigma=True flag. Unten ist die Ausgabe, die

1. diagonale Elemente der Kovarianzmatrix cov (0,0) und
2. cov (1,1),
3. falscher Weg für Standardabweichung von den Ausgängen für Steigung und
4. falscher Weg für die Konstante und
5. rechts Weg für Standardabweichung von den Ausgängen für Steigung und
6. rechts Weg für die Konstante

odr: 1.739631e-06 0.02302262 [ 0.00014863 0.0170987 ] [ 0.00131895 0.15173207] curve_fit: 2.209469e-08 0.00029239 [ 0.00014864 0.01709943] [ 0.0004899 0.05635713] polyfit: 2.232016e-08 0.00029537 [ 0.0001494 0.01718643]

Beachten Sie, dass odr und polyfit genau die gleiche Standardabweichung haben. Polyfit berechnet die Unsicherheiten als Eingabe, so dass odr keine Unsicherheiten bei der Berechnung der Standardabweichung verwendet . Die Kovarianzmatrix verwendet diese und wenn im odr-Fall die Standardabweichung aus der Kovarianzmatrix berechnet wird, gibt es Unsicherheiten, die sich ändern, wenn die Unsicherheit erhöht wird. Fummeln mit dy im folgenden Code wird es zeigen.

Ich schreibe das hier vor allem, weil dies wichtig ist, wenn man Fehlergrenzen herausfindet (und der fortran odrpack guide, auf den scipy verweist, enthält irreführende Informationen: Die Standardabweichung sollte die Quadratwurzel der Kovarianzmatrix sein, wie es in der Anleitung steht aber es ist nicht).