Python / Clé basse automatique inadaptée aux données expérimentales
1.Tout d'abord
Je l'ai écrit dans le but d'ajuster un grand nombre de données expérimentales à la fois. Il est difficile d'ouvrir les données une par une, de les ajuster, d'enregistrer les paramètres du résultat, etc. Prenons Lawrench Anne comme exemple, mais je pense que le flux lui-même de création d'une liste de fichiers → ajustement peut avoir d'autres applications.
2. Histoire facile de Lawrench Anne (skip OK)
Lorenzian (fonction de Lorentz) est une telle fonction.
3. Exemples de données
Je pense qu'il est plus facile de comprendre si vous jouez réellement avec les données, je vais donc créer des exemples de données.
dataset.py
import numpy as np
#Définition des paramètres
intensity = 50 #Force
HWHM = 5 #Demi-prix demi-largeur
a = 10 #Grande quantité de variation de données
#Créer un fichier de données
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')
Lors de l'exécution, 7 fichiers dat seront créés dans le répertoire, séparés par 10 de X0 = -30 à 30. Données numériques sur deux colonnes (X, Y) séparées par des tabulations. X, Y (chaîne de caractères) est écrit comme un index sur la première ligne. Le contenu est la fonction Lorentz de plus tôt, Y = $ f ($ X $) $. Y reçoit une légère variation avec un nombre aléatoire. A titre d'exemple, le tracé de "X0 = 0,00at" ressemble à celui ci-dessous.
À partir de maintenant, je voudrais mettre à mal ces 7 fichiers dat à la fois pour créer un fichier dat qui résume le graphique de chaque résultat et les paramètres convergés.
4. En vrac inadapté
Le sujet principal est ci-dessous. Ajustez les 7 fichiers dat créés à la fois et n'obtenons que le résultat. Je voudrais utiliser jupyter pour que la procédure soit facile à suivre. Si vous écrivez d'abord le flux général à partir d'ici,
- Créez une liste de fichiers adaptés
- Ajustement des fichiers dans la liste un par un avec itération
- Dessinez le graphique et écrivez les paramètres d'ajustement convergés dans un autre fichier
Sera.
4.1 Faire une liste de fichiers pour s'adapter
Utilisez le module glob pour créer une liste de fichiers à ajuster, filelist, pour traitement par l'itérateur. Ici, le caractère générique «*» est utilisé pour lister uniquement les fichiers nommés «X0 = ○○ .dat» dans le dossier. Alors je vais le faire avec Jupyter,
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']
Référence: Python: Trier par numéro de nom de fichier (Site externe)
Éxécuter. Ce que fait ʻIn [2] est de créer une liste et de trier les fichiers dans la liste. Si vous l'attrapez simplement avec glob`, il sera trié par chaîne de caractères, donc il ne sera pas trié par ordre croissant de nombres. Ce n'est pas un gros problème, mais je me sens toujours mal à l'aise plus tard, donc je les trie par ordre croissant.
4.2 Estimation et ajustement des paramètres initiaux
L'intérieur de la boucle devient un peu plus long, mais cela ressemble à ceci. (J'expliquerai plus tard.)
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]))
#Estimation de la valeur initiale (intensité et valeur centrale)
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)]
#Réglage de la valeur initiale
pini = np.array([intensity, X0, 1])
#raccord
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))
#Voir les résultats
print("initial parameter\toptimized parameter")
for i, v in enumerate(pini):
print(str(v)+ '\t' + str(popt[i]) + ' ± ' + str(perr[i]))
#Écriture du résultat de l'ajustement dans le fichier dat
result.writelines(name + '\t' + '\t'.join([str(p) + '\t' + str(e) for p, e in zip(popt, perr)]) + '\n')
#Créer un tableau pour dessiner une courbe d'ajustement
fitline = func(X, popt[0], popt[1], popt[2])
#Tracé des résultats et sauvegarde de l'image
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]Résultat de
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
・ ・ ・
(Omis)
・ ・ ・
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
Je voudrais expliquer chacun d'eux. Tout d'abord, bien qu'il soit supérieur à l '# estimation de la valeur initiale commentée, nous préparons un fichier pour écrire le résultat final (création d'un fichier et écriture d'un index). Ensuite, l'itération est lancée à l'aide de la liste des fichiers de données sous l'instruction for.
python
match = re.match('(.*).dat', filename)
name = match.group(1)
À la place, la chaîne de caractères avant l'extension de .dat est extraite avec nom pour le nom de fichier lors de l'exportation de l'image à la fin.
Vient ensuite l '«# estimation de la valeur initiale». Dans cet ajustement, les ** paramètres initiaux ** (pini) doivent être ** intensité de crête ** (ʻintensité), ** X ** ( X0) lors de la prise d'un pic, ** moitié prix. Il est demi-largeur ** (HWHM`). Il est difficile d'estimer avec quelle précision, mais ici nous estimons chacun comme suit.
Dans l'ordre d'estimation,
ʻIntensity: extrait les valeurs maximale et minimale de Y et adopte celle avec la valeur absolue la plus élevée. X0: adopte la valeur de X lorsque l''intensité adoptée apparaît
«HWHM»: non estimé. Je l'ai mis à 1 pour le moment. Les deux paramètres ci-dessus sont assez précis, il semble donc bon de faire un compromis ici.
Et au «# essayage». Ajuster en utilisant les paramètres initiaux estimés. Ici, le montage se fait selon la procédure de montage de SciPy. Les paramètres convergents sont stockés dans popt. perr est l'erreur standard. Il prend la racine de la composante diagonale de la covariance «pcov» des paramètres d'ajustement.
référence:
- Modélisation de fonctions non linéaires avec Python ・ Scipy.optimize.curve_fit (document officiel SciPy)
Le dernier est «# Afficher le résultat» et «# Écrire le résultat de l'ajustement dans le fichier dat», «# Créer un tableau pour dessiner la courbe d'ajustement», «# Tracer le résultat et enregistrer l'image». Les # affichage du résultat et # écriture du résultat ajusté dans le fichier dat sont omis tels quels. «#Créer un tableau pour dessiner la courbe d'ajustement» est le résultat de la substitution de X pour les données dans la fonction définie à l'aide des paramètres convergés. Donc, ici, les données de la ligne d'ajustement sont préparées avec le même score que les données d'origine. Le # tracé de résultat et image de sauvegarde est également omis tel quel.
python
plt.savefig(name + ".png ")
De cette façon, j'utilise le "nom" qui est le résultat du premier "re.match".
5. Résultat
L'image du résultat approprié ressemble à ceci.
Ensuite, le résultat numérique (les paramètres d'ajustement convergés et l'erreur standard sont délimités par des tabulations) est écrit dans fitresult.dat.
6. Conclusion
Ainsi, j'ai pu adapter 7 fichiers à la fois. Cette fois, il était sept pour des raisons de simplicité, mais si vous faites une expérience, vous devrez peut-être adapter des dizaines à plus d'une centaine de données dans certains cas, donc je pense que ce sera utile dans de tels cas.
7. Supplément
Si le pic est si clair, n'est-il pas possible de fixer correctement tous les paramètres initiaux? Je pense que certaines personnes le pensent. Cependant, s'il y a un écart dans les paramètres initiaux au niveau de la commande,
La réalité est que ce sera quelque chose comme ça. Je ne pense pas qu'il soit nécessaire d'estimer la valeur initiale de l'ajustement linéaire, mais si la fonction devient un peu compliquée, il semble préférable d'estimer la valeur initiale dans une certaine mesure (car il suffit d'estimer l'ordre).
Certaines personnes peuvent ne pas être familières avec Jupyter, je vais donc mettre celui sous la forme de .py à la fin.
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 = []
'''
/////////////////////////////////
Tri des fichiers de données
/////////////////////////////////
'''
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 (clé basse inadaptée)
/////////////////////////////////
'''
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]))
#Estimation de la valeur initiale (intensité et valeur centrale)
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)]
#Réglage de la valeur initiale
pini = np.array([intensity, X0, 1])
#raccord
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))
#Voir les résultats
print("initial parameter\toptimized parameter")
for i, v in enumerate(pini):
print(str(v)+ '\t' + str(popt[i]) + ' ± ' + str(perr[i]))
#Écriture du résultat de l'ajustement dans le fichier dat
result.writelines(name + '\t' + '\t'.join([str(p) + '\t' + str(e) for p, e in zip(popt, perr)]) + '\n')
#Créer un tableau pour dessiner une courbe d'ajustement
fitline = []
for v in X:
fitline.append(func(v, popt[0], popt[1], popt[2]))
#Tracé des résultats et sauvegarde de l'image
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