[PYTHON] Vergleich der Anpassungsprogramme

Einführung

Um ein Handbuch für das Analyseverfahren und ein Beispielprogramm für das Experiment eines Studenten im Grundstudium zu schreiben, habe ich versucht, mithilfe verschiedener Programme herauszufinden, wie geeignete Daten mit der Minimum-Square-Methode angepasst werden können, und werde daher eine Notiz davon hinterlassen.

Da das Prinzip insofern dasselbe ist, als es die Methode der kleinsten Quadrate verwendet, ist es ein Problem, wenn Sie nicht das gleiche Ergebnis erhalten, aber leider wusste ich vorher, dass die Antwort je nach Programm unterschiedlich ist, aber in letzter Zeit die Ursache und das Mittel Ich habe den Text gefunden, der besagt, also habe ich beschlossen, ihn für mich selbst zu organisieren. Dies ist der Text, der es ausgelöst hat. Peter Young, "Everything you wanted to know about Data Analysis and Fitting but were afraid to ask" Anscheinend scheint es sich um ein Vorlesungsmaterial einer Universität zu handeln. Hier sagen Gnuplot und Pythons `` `scipy.optimize.curve_fit```, dass bei einem Fehler in den Daten der an den resultierenden Parameter angehängte Fehlerwert falsch ist und korrigiert werden muss.

Im Fall von Gnuplot ist es notwendig, es so zu ändern.

to get correct error bars on fit parameters from gnuplot when there are error bars on the points, you have to divide gnuplot’s asymptotic standard errors by the square root of the chi-squared per degree of freedom (which gnuplot calls FIT STDFIT and, fortunately, computes correctly).

Die Python `scipy.optimize``` ist kompliziert, die` kurve_fit muss geändert werden, aber die `` `Leastsq muss nicht geändert werden.

I recently learned that error bars on fit parameters given by the routine curve_fit of python also have to be corrected in the same way. This is shown in two of the python scripts in appendix H. Curiously, a different python fitting routine, leastsq, gives the error bars correctly.

Es heißt, ich frage mich, warum sich niemand dessen bewusst ist.

It is curious that I found no hits on this topic when Googling the internet.

Beispieldaten

Ich habe einen solchen Datensatz vorbereitet. 20 Sätze (x, y, ey). Ich werde dies in einer geraden Linie passen. Ich werde es als Textdatei `` `data.txt``` speichern.

0.0e+00 -4.569987720595017344e-01 1.526828747143463172e+00
1.0e+00 -7.106162255269843353e-01 1.402885069270964458e+00
2.0e+00 1.105159634902675325e+00 1.735638554786020915e+00
3.0e+00 -1.939878950652441869e+00 1.011014634823069747e+00
4.0e+00 3.609690931525689983e+00 1.139915698020605550e+00
5.0e+00 8.535035219721383015e-01 9.338187791237286817e-01
6.0e+00 4.770810591544029755e+00 1.321364026236713451e+00
7.0e+00 3.323982457761388787e+00 1.703973901689593173e+00
8.0e+00 3.100622722027332578e+00 1.002313080286136637e+00
9.0e+00 4.527766245564444070e+00 9.876090792441625243e-01
1.0e+01 1.990062497396323682e+00 1.355607177365929505e+00
1.1e+01 5.113013340421659336e+00 9.283045349565146598e-01
1.2e+01 4.391676777018354905e+00 1.337677147217683160e+00
1.3e+01 5.388022504497612886e+00 9.392443558621643707e-01
1.4e+01 1.134921361159764075e+01 9.232583484294124565e-01
1.5e+01 6.067025020573844074e+00 1.186258237028150475e+00
1.6e+01 1.052771612360148445e+01 1.200732350014090954e+00
1.7e+01 6.221953870216905713e+00 8.454085761899273743e-01
1.8e+01 9.628358150028700990e+00 1.442970173161927772e+00
1.9e+01 9.493784288063746857e+00 8.196526623903285236e-01

Beispielcode in verschiedenen Programmen

CERN ROOT

Das ist ROOT. https://root.cern.ch Der Code war der kürzeste. Im Fall von ROOT wird es, wenn es einfach ist, vorbereitet, ohne eine Anpassungsfunktion zu definieren.

fit.C

{
  TGraphErrors *g = new TGraphErrors("data.txt","%lg %lg %lg");
  g->Fit("pol1");
  g->Draw("ape");
}

Das Ergebnis sieht so aus.

Gnuplot

Gnuplot passt auch gut zu Funktionen, sodass der erforderliche Code sehr kurz ist. Ich denke, das Problem mit Gnuplot ist, dass die Standardausgabe zu schlecht aussieht.

fit.gp

set fit errorvariables

f(x) = p0 + p1*x
fit f(x) "data.txt" u 1:2:3 via p0,p1

plot "data.txt" u 1:2:3 w yerr, f(x)

print "\n ====== ERROR CORRECTED ========"
print "Chi2/NDF = ",FIT_STDFIT**2 * FIT_NDF,"/",FIT_NDF
print "  p0 = ",p0," +- ",p0_err/FIT_STDFIT
print "  p1 = ",p1," +- ",p1_err/FIT_STDFIT

# ---------Das Folgende ist die Aussehensanpassung der Figur--------

set term qt font "Helvetica"
set xrange [-1:20]
set rmargin 5
set tics font ",16"
set key font ",16"
set key left top
set bars fullwidth 0
set style line 1 lc "#0080FF" lw 1 pt 7 ps 1
set style line 2 lc "#FF3333" lw 2 pt 0 ps 1

set label 1 at first 1,11 sprintf("Chi2/ndf = %5.2f / %2d",FIT_STDFIT**2 * FIT_NDF,FIT_NDF) font ",18"
set label 2 at first 1,10 sprintf("p0 = %6.3f +- %7.4f",p0,p0_err/FIT_STDFIT) font ",18"
set label 3 at first 1,9  sprintf("p1 = %6.4f +- %7.5f",p1,p1_err/FIT_STDFIT) font ",18"

plot "data.txt" u 1:2:3 w yerr ls 1,\
     f(x) ls 2

Der in der Standardausgabe geschriebene "Asymptotic Standard Error" ist falsch und muss korrigiert werden. Teilen Sie den Fehlerwert insbesondere durch eine Variable namens FIT_STDFIT, wie im obigen Code. Wenn Sie zu Beginn `` `set fit errorvariables``` schreiben, können Sie den Fehlerwert auch mit dem Variablennamen _err abrufen. Wenn Sie es ändern, wird derselbe Wert wie ROOT angezeigt.

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
p0              = -1.06859         +/- 0.9578       (89.64%)
p1              = 0.566268         +/- 0.07983      (14.1%)

correlation matrix of the fit parameters:
                p0     p1     
p0              1.000 
p1             -0.884  1.000 

 ====== ERROR CORRECTED ========
Chi2/NDF = 59.1533703771407/18
  p0 = -1.06858871709936 +- 0.528376469987239
  p1 = 0.566267669300731 +- 0.0440357299923021

Beachten Sie, dass Sie diesen Fehler nur beheben sollten, wenn ein Fehler im Datenpunkt vorliegt (wenn Sie beim Aufnehmen von Daten mit dem Befehl fit drei Spalten nach `` `using``` angeben). Was musst du machen. Nehmen Sie diese Korrektur nicht vor, wenn alle mit dem gleichen Gewicht passen (= Daten ohne Fehler).

python scipy.optimize.curve_fit

Das Laden der Daten war mit `` `numpy.loadtxt``` einfach.

fit_curve_fit.py

#Daten lesen
import numpy as np

data = np.loadtxt("data.txt")
xx = data.T[0] 
yy = data.T[1]
ey = data.T[2]

#Definieren Sie eine Anpassungsfunktion ff
def ff(x,a,b):
    return a + b*x

#Ergebnisse anpassen und anzeigen
from scipy.optimized import curve_fit
import math

par, cov = curve_fit(ff,xx,yy,sigma=ey)

chi2 = np.sum(((func_pol1(xx,par[0],par[1])-yy)/ey)**2)
print("chi2 = {:7.3f}".format(chi2))
print("p0 : {:10.5f} +- {:10.5f}".format(par[0],math.sqrt(cov[0,0]/chi2*18)))
print("p1 : {:10.5f} +- {:10.5f}".format(par[1],math.sqrt(cov[1,1]/chi2*18)))

#Anzeige in der Grafik
import matplotlib.pyplot as plt

x_func = np.arange(0,20,0.1)
y_func = par[0] + par[1]*x_func

plt.errorbar(xx,yy,ey,fmt="o")
plt.plot(x_func,y_func)
plt.show()

Es scheint, dass Gnuplots FIT_STDFIT nicht bereitgestellt wird, daher werde ich Chi2 und NDF selbst berechnen und den Parameterfehler unter Verwendung der diagonalen Komponente der Ausgabekovarianzmatrix berechnen. Wenn Sie es richtig berechnen, erhalten Sie den richtigen Wert.

chi2 =  59.153
p0 :   -1.06859 +-    0.52838
p1 :    0.56627 +-    0.04404

python scipy.optimize.leastsq

Ich habe das also noch nie benutzt https://qiita.com/yamadasuzaku/items/6d42198793651b91a1bc Ich durfte mich beziehen. Es war ein wenig verwirrend (nicht Chi ^ 2), dass ich Chi vorbereiten musste, nicht die Funktion, die ich anpassen wollte.

fit_leastsq.py

#Daten lesen
import numpy as np

data = np.loadtxt("data.txt")
xx = data.T[0] 
yy = data.T[1]
ey = data.T[2]

#Definiere Chi
from scipy.optimize import leastsq
import math

def chi(prm,x,y,ey):
    return (((prm[0]+prm[1]*x)-y)/ey)

#Bereiten Sie den Anfangswert vor und passen Sie ihn an
init_val = (-0.5, 0.5)

prm, cov, info, msg, ier = leastsq(chi,init_val,args=(xx,yy,ey),full_output=True)

chi2 = np.sum((((prm[0]+prm[1]*xx) - yy)/ey)**2)

print("chi2 = {:7.3f}".format(chi2))
print("p0 : {:10.5f} +- {:10.5f}".format(prm[0],math.sqrt(cov[0,0])))
print("p1 : {:10.5f} +- {:10.5f}".format(prm[1],math.sqrt(cov[1,1])))

Die Grafikanzeige ist dieselbe wie oben, daher wird sie weggelassen. Das Ergebnis ist wie folgt. Im Fall von Leastsq ist keine Modifikation erforderlich, sodass die Quadratwurzel der Diagonalkomponente der Ausgangskovarianzmatrix unverändert verwendet werden kann.

chi2 =  59.153
p0 :   -1.06859 +-    0.52838
p1 :    0.56627 +-    0.04404

schließlich

Früher dachte ich, dass das Ergebnis von Gnuplot möglicherweise nicht mit der manuellen Berechnung übereinstimmt, aber ich verwende normalerweise ROOT, also habe ich es nicht ernsthaft überprüft, aber ich habe schließlich herausgefunden, wie ich damit umgehen soll, also war es sehr erfrischend. Persönlich bin ich daran gewöhnt, daher ist ROOT einfach, aber ich dachte, dass Gnuplot oder Pythons Curve_fit leichter zu verstehen wären, wenn ich Studenten unterrichten würde. Beide haben jedoch das Problem, dass der Fehler, der an den resultierenden Parameter angehängt werden sollte, korrigiert werden muss.

Übrigens dachte ich, dass es heutzutage besser wäre, Python anstelle von C oder Gnuplot für Studenten zu unterrichten, also dachte ich, dass es besser wäre, Python und mein eigenes Studium zu unterrichten. Ich habe es versucht. Sicher, Python ist insofern gut, als es alles von der Datenverarbeitung bis zur Diagrammanzeige kann, aber wenn es um Funktionsdefinition und Diagrammanzeige geht, ist Gnuplot intuitiver und auf das Zeichnen von Diagrammen spezialisiert. Ich hatte auch das Gefühl, dass es nur einen gab. Als einfaches Beispiel zeigt es fast dasselbe an, aber wenn Sie die beiden unten vergleichen, denke ich, dass Gnuplot intuitiver ist. Python sieht allerdings besser aus.

gnuplot

set xrange [0:10]
f(x) = sin(x)
plot f(x)

python

import numpy as np
import matplotlib.pyplot as plt
x = np.arange(0,10,0.1)
y = np.sin(x)
plt.plot(x,y)
plt.show()