Python / Automatic Low Wrench, nicht für experimentelle Daten geeignet
1. Zuallererst
Ich habe es geschrieben, um viele experimentelle Daten gleichzeitig anzupassen. Es ist mühsam, die Daten einzeln zu öffnen, anzupassen, die Parameter des Ergebnisses aufzuzeichnen und so weiter. Nehmen wir als Beispiel Lawrench Anne, aber ich denke, dass der Ablauf des Erstellens einer Dateiliste → Anpassung andere Anwendungen haben kann.
2. Einfache Geschichte von Lawrench Anne (OK überspringen)
Lorentian (Lorentz-Funktion) ist eine solche Funktion.
3. Beispieldaten
Ich denke, es ist einfacher zu verstehen, ob Sie tatsächlich mit den Daten spielen, also werde ich Beispieldaten erstellen.
dataset.py
import numpy as np
#Parameterdefinition
intensity = 50 #Stärke
HWHM = 5 #Halber Preis, halbe Breite
a = 10 #Große Menge an Datenvariationen
#Erstellen einer Datendatei
for X0 in range (-30, 31, 10):
filename = 'X0=' + str(X0) + '.dat'
with open(filename, 'w') as file:
file.writelines('X' + '\t' + 'Y' +'\n')
for X in range (-100,101):
Y = intensity * HWHM**2 / ((X-X0)**2 + HWHM**2) + a * np.random.rand() - 5
file.writelines(str(X) + '\t' + str(Y) + '\n')
Bei der Ausführung werden 7 Datendateien im Verzeichnis erstellt, die durch 10 von X0 = -30 bis 30 getrennt sind. Numerische Daten in zwei durch Tabulatoren getrennten Spalten (X, Y). X, Y (Zeichenkette) wird als Index in die erste Zeile geschrieben. Der Inhalt ist die Lorentz-Funktion von früher, Y = $ f ($ X $) $. Y erhält eine leichte Variation mit einer Zufallszahl. Das Zeichnen von "X0 = 0.00at" sieht beispielsweise wie folgt aus.
Von nun an möchte ich diese 7 Datendateien auf einmal nicht mehr anpassen, um eine Datendatei zu erstellen, die das Diagramm jedes Ergebnisses und die konvergierten Parameter zusammenfasst.
4. Bulk nicht passend
Das Hauptthema ist unten. Passen wir die erstellten 7 Datendateien gleichzeitig an und erhalten nur das Ergebnis. Ich möchte jupyter verwenden, damit das Verfahren leicht zu befolgen ist. Wenn Sie zuerst den allgemeinen Ablauf von hier aus schreiben,
- Erstellen Sie eine Liste der passenden Dateien
- Anpassen der Dateien in der Liste nacheinander mit Iteration
- Zeichnen Sie das Diagramm und schreiben Sie die konvergierten Anpassungsparameter in eine andere Datei
Wird sein.
4.1 Erstellen Sie eine Liste der passenden Dateien
Verwenden Sie das Glob-Modul, um eine Liste von Dateien zu erstellen, die zur Verarbeitung durch den Iterator passen (Dateiliste). Hier wird der Platzhalter "*" verwendet, um nur die Dateien mit dem Namen "X0 = ○○ .dat" im Ordner aufzulisten. Also werde ich es mit jupyter machen,
In[1]
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ptick
import scipy.optimize
import glob, re
from functools import cmp_to_key
In[2]
filelist = glob.glob('X0=*.dat')
filelist2 = []
for filename in filelist:
match = re.match('X0=(.*).dat', filename)
tuple = (match.group(1), filename)
filelist2.append(tuple)
def cmp(a, b):
if a == b: return 0
return -1 if a < b else 1
def cmptuple(a, b):
return cmp(float(a[0]), float(b[0]))
filelist = [filename[1] for filename in sorted(filelist2, key=cmp_to_key(cmptuple))]
filelist
Out[2]
['X0=-30.dat',
'X0=-20.dat',
'X0=-10.dat',
'X0=0.dat',
'X0=10.dat',
'X0=20.dat',
'X0=30.dat']
Referenz: Python: Nach Dateinamen sortieren (Externe Site)
Ausführen. In [2] erstellen Sie eine Liste und sortieren die Dateien in der Liste. Wenn Sie es nur mit "glob" abfangen, wird es nach Zeichenfolgen sortiert, sodass es nicht in aufsteigender Reihenfolge der Zahlen sortiert wird. Das ist kein großes Problem, aber ich fühle mich später immer noch unwohl, also sortiere ich sie in aufsteigender Reihenfolge.
4.2 Anfängliche Parameterschätzung und -anpassung
Das Innere der Schleife wird auf jeden Fall etwas länger, aber es sieht so aus. (Ich werde es später erklären.)
In[3]
result = open('fitting_result.dat', 'w')
result.writelines("filename\tA\tA_err\tX0\tX0_err\tHWHM\tHWHM_err\n")
for filename in filelist:
print(filename)
match = re.match('(.*).dat', filename)
name = match.group(1)
with open(filename, 'r') as file:
X, Y = [], []
for line in file.readlines()[1:]:
items = line[:-1].split('\t')
X.append(float(items[0]))
Y.append(float(items[1]))
#Anfangswertschätzung (Intensität und Mittelwert)
max_Y = max(Y)
min_Y = min(Y)
if np.abs(max_Y) >= np.abs(min_Y):
intensity = max_Y
X0 = X[Y.index(max_Y)]
else:
intensity = min_Y
X0 = X[Y.index(min_Y)]
#Anfangswerteinstellung
pini = np.array([intensity, X0, 1])
#passend zu
def func(x, intensity, X0, HWHM):
return intensity * HWHM**2 / ((x-X0)**2 + HWHM**2)
popt, pcov = scipy.optimize.curve_fit(func, X, Y, p0=pini)
perr = np.sqrt(np.diag(pcov))
#Ergebnisse anzeigen
print("initial parameter\toptimized parameter")
for i, v in enumerate(pini):
print(str(v)+ '\t' + str(popt[i]) + ' ± ' + str(perr[i]))
#Schreiben des Anpassungsergebnisses in die Datendatei
result.writelines(name + '\t' + '\t'.join([str(p) + '\t' + str(e) for p, e in zip(popt, perr)]) + '\n')
#Erstellen Sie ein Array zum Zeichnen einer Anpassungskurve
fitline = func(X, popt[0], popt[1], popt[2])
#Ergebnisplot und Bild speichern
plt.clf()
plt.figure(figsize=(10,6))
plt.plot(X, Y, 'o', alpha=0.5)
plt.plot(X, fitline, 'r-', linewidth=2)
plt.xlabel('X')
plt.ylabel('Y')
plt.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
plt.savefig(name + ".png ")
plt.close()
result.close()
In[3]Ergebnis von
X0=-30.dat
initial parameter optimized parameter
52.8746388447 52.1363835425 ± 1.37692592137
-29.0 -30.2263312447 ± 0.132482308868
1.0 5.01691953143 ± 0.187432386079
・ ・ ・
(Weggelassen)
・ ・ ・
X0=30.dat
initial parameter optimized parameter
50.0825336071 50.5713312616 ± 1.54634448135
31.0 30.1170389209 ± 0.149532401478
1.0 4.89078858649 ± 0.211548017933
Ich möchte jeden einzelnen erklären. Erstens bereiten wir, obwohl es über der auskommentierten "# Schätzung des Anfangswertes" liegt, eine Datei vor (Erstellen einer Datei und Schreiben eines Index), um zuerst das Endergebnis zu schreiben. Anschließend wird die Iteration mithilfe der Liste der Datendateien unterhalb der for-Anweisung gestartet.
python
match = re.match('(.*).dat', filename)
name = match.group(1)
An der Stelle wird die Zeichenfolge vor der Erweiterung von ".dat" mit "name" für den Dateinamen extrahiert, wenn das Bild am Ende exportiert wird.
Als nächstes folgt die "# Schätzung des Anfangswertes". Bei dieser Anpassung müssen ** Anfangsparameter ** (pini) ** Spitzenintensität ** ( Intensität), ** X ** (X0) sein, wenn ein Peak genommen wird, ** halber Preis. Es ist halb so breit ** (HWHM). Es ist schwierig zu schätzen, mit welcher Genauigkeit, aber hier schätzen wir jede wie folgt.
In der Reihenfolge der Schätzung,
Intensität: Extrahiert die Maximal- und Minimalwerte von Y und übernimmt den mit dem größeren Absolutwert.
"X0": Übernehmen Sie den Wert von X, wenn die angenommene "Intensität" angezeigt wird
HWHM: Nicht geschätzt. Ich habe es vorerst auf 1 gesetzt. Die beiden oben genannten Parameter sind ziemlich genau, daher scheint es gut, hier Kompromisse einzugehen.
Und bei der "# Anpassung". Passen Sie mit den geschätzten Anfangsparametern an. Hier erfolgt die Anpassung nach dem Anpassungsverfahren von SciPy. Konvergierte Parameter werden in "popt" gespeichert. perr ist der Standardfehler. Es ist die Wurzel der diagonalen Komponente der Kovarianz "pcov" der Anpassungsparameter.
Referenz:
- Nichtlineare Funktionsmodellierung mit Python ・ Scipy.optimize.curve_fit (offizielles SciPy-Dokument)
Das letzte ist "# Ergebnis anzeigen" und "Anpassungsergebnis in Datendatei schreiben", "Array zum Zeichnen der Anpassungskurve erstellen", "Ergebnis zeichnen und Bild speichern". Die "# Anzeige des Ergebnisses" und "# Schreiben des Anpassungsergebnisses in die Datendatei" werden so wie sie sind weggelassen. #Erstellen eines Arrays zum Zeichnen der Anpassungskurve ist das Ergebnis des Ersetzens der Daten durch X in die Funktion, die mithilfe der konvergierten Parameter definiert wurde. Daher werden hier die Anpassungsliniendaten mit der gleichen Punktzahl wie die Originaldaten erstellt. Das # Ergebnisdiagramm und Bildspeicherung wird ebenfalls weggelassen.
python
plt.savefig(name + ".png ")
Auf diese Weise verwende ich den "Namen", der das Ergebnis des ersten "re.match" ist.
5. Ergebnis
Das Bild des Anpassungsergebnisses sieht so aus.
Dann wird das numerische Ergebnis (konvergierte Anpassungsparameter und Standardfehler werden durch Tabulatoren getrennt) in "fitresult.dat" geschrieben.
6. Fazit
So konnte ich 7 Dateien gleichzeitig anpassen. Diesmal war es der Einfachheit halber sieben, aber wenn Sie ein Experiment durchführen, müssen Sie möglicherweise Dutzende an mehr als hundert Daten anpassen, sodass ich denke, dass dies in solchen Fällen nützlich sein wird.
7. Ergänzung
Wenn der Peak so klar ist, ist es nicht möglich, alle Anfangsparameter entsprechend festzulegen? Ich denke, einige Leute denken das. Wenn es jedoch zu einer Abweichung der Anfangsparameter auf Auftragsebene kommt,
Die Realität ist, dass es so etwas sein wird. Ich denke nicht, dass es notwendig ist, den Anfangswert für die lineare Anpassung zu schätzen, aber wenn die Funktion etwas kompliziert wird, scheint es besser, den Anfangswert bis zu einem gewissen Grad zu schätzen (da es ausreicht, die Reihenfolge zu schätzen).
Einige Leute kennen Jupyter vielleicht nicht, deshalb werde ich die am Ende in Form von ".py" setzen.
lorentzian_fittng.py
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
import glob, re
from functools import cmp_to_key
filelist = glob.glob('X0=*.dat')
filelist2 = []
'''
/////////////////////////////////
Datendateien sortieren
/////////////////////////////////
'''
for filename in filelist:
match = re.match('X0=(.*).dat', filename)
tuple = (match.group(1), filename)
filelist2.append(tuple)
def cmp(a, b):
if a == b: return 0
return -1 if a < b else 1
def cmptuple(a, b):
return cmp(float(a[0]), float(b[0]))
filelist = [filename[1] for filename in sorted(filelist2, key=cmp_to_key(cmptuple))]
'''
/////////////////////////////////
Analyse (niedriger Schraubenschlüssel nicht passend)
/////////////////////////////////
'''
result = open('fitting_result.dat', 'w')
result.writelines("filename\tA\tA_err\tX0\tX0_err\tHWHM\tHWHM_err\n")
for filename in filelist:
print(filename)
match = re.match('(.*).dat', filename)
name = match.group(1)
with open(filename, 'r') as file:
X, Y = [], []
for line in file.readlines()[1:]:
items = line[:-1].split('\t')
X.append(float(items[0]))
Y.append(float(items[1]))
#Anfangswertschätzung (Intensität und Mittelwert)
max_Y = max(Y)
min_Y = min(Y)
if np.abs(max_Y) >= np.abs(min_Y):
intensity = max_Y
X0 = X[Y.index(max_Y)]
else:
intensity = min_Y
X0 = X[Y.index(min_Y)]
#Anfangswerteinstellung
pini = np.array([intensity, X0, 1])
#passend zu
def func(x, intensity, X0, HWHM):
return intensity * HWHM**2 / ((x-X0)**2 + HWHM**2)
popt, pcov = scipy.optimize.curve_fit(func, X, Y, p0=pini)
perr = np.sqrt(np.diag(pcov))
#Ergebnisse anzeigen
print("initial parameter\toptimized parameter")
for i, v in enumerate(pini):
print(str(v)+ '\t' + str(popt[i]) + ' ± ' + str(perr[i]))
#Schreiben des Anpassungsergebnisses in die Datendatei
result.writelines(name + '\t' + '\t'.join([str(p) + '\t' + str(e) for p, e in zip(popt, perr)]) + '\n')
#Erstellen Sie ein Array zum Zeichnen einer Anpassungskurve
fitline = []
for v in X:
fitline.append(func(v, popt[0], popt[1], popt[2]))
#Ergebnisplot und Bild speichern
plt.clf()
plt.figure(figsize=(10,6))
plt.plot(X, Y, 'o', alpha=0.5)
plt.plot(X, fitline, 'r-', linewidth=2)
plt.xlabel('X')
plt.ylabel('Y')
plt.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
plt.savefig(name + ".png ")
plt.close()
result.close()
Recommended Posts