[PYTHON] Analyse des données de mesure (2) -Hydrobacter et raccord, recommandation lmfit-

Préface

Rencontre d'un jour, Kyo "T-kun, ne peux-tu pas ajuster un peu mieux ce graphique?" "Ne pouvons-nous pas changer la gamme d'ajustement ou la pondérer?" T "... je vais essayer"

Il veut modifier le raccord. Il est souvent nécessaire de changer la procédure d'analyse sur un coup de tête, et il est difficile d'y faire face. Je veux pouvoir le faire sur place.

Cette fois, au lieu de scipy.optimize.curve_fit, je vais introduire une nouvelle bibliothèque appelée lmfit. L'ajustement à l'aide de scipy a été effectué dans Analyse des données de mesure ①-Mémorandum de raccord scipy-. Le codage est fait sur Jupyter & Jupyter recommandé.

J'espère que je peux bien présenter lmfit.

GithHub dispose également d'un bloc-notes et d'exemples de données. Depuis ici

À propos de lmfit

lmfit est développé avec une compatibilité ascendante avec scipy.optimize. Vous pouvez contraindre les paramètres et les erreurs de pondération lors de l'ajustement. Il semble que cela fonctionne comme ajouter diverses informations à l'objet et enfin exécuter la commande fit pour l'optimiser. Si vous vous y habituez, ce sera plus pratique que scipy.

Cette expérience

Les données d'impulsion ont été obtenues à partir de l'instrument de mesure. Commencez par créer un histogramme avec la valeur de crête de l'impulsion. Après cela, l'ajustement est effectué et la netteté du spectre est évaluée. De plus, nous essayerons d'améliorer la netteté en utilisant l'option de pondération lmfit.

Lire les données

Lisez les données d'impulsion balayées de l'instrument de mesure avec des pandas. Encore une fois, utilisez tkinter pour connecter les chemins de fichiers.

.ipynb

import tkinter 
from tkinter import filedialog
import numpy as np
import pandas as pd

root = tkinter.Tk()
root.withdraw()
file = filedialog.askopenfilename(
title='file path_',
filetypes=[("csv","csv")])
if file:
    pulse = pd.read_csv(file,engine='python',header=0,index_col=0)

pulse.iloc[:,[0,1,2]].plot()

Tracons seulement 3 données. Il a été balayé de sorte que la valeur de crête maximale puisse être obtenue près du 50e point.

Créer un histogramme

Créez un histogramme avec la valeur de crête maximale. Lors de la création d'un histogramme, vous devez déterminer la taille du pas sur l'axe horizontal. Cette fois, selon la formule de Sturges, le nombre de cases est fixé à: k = 1 + Log2 (n).

* Np.histogram renvoie ** width ** autant que le nombre.

.ipynb

bins = int(1+np.log2(pulse.shape[1]))#Officiel de Sturges
count, level = np.histogram(pulse.max(),bins=bins)
level = np.array(list(map(lambda x:(level[x]+level[x+1])/2,range(level.size-1))))
#Aligner la fréquence et la taille
level = np.array(list(map(lambda x:(level[x]+level[x+1])/2,range(level.size-1))))
import matplotlib.pyplot as plt
plt.plot(level,count,marker='.',linestyle='None')
plt.show()

Il y avait une petite montagne sur le côté droit. Séparons les données après 50 et effectuons à nouveau une distribution pour la montagne sur le côté droit.

.ipynb

mask = pulse.max()>50#Créer un tableau de valeurs booléennes
pulse = pulse.loc[:,mask]#Masquage des données d'impulsion

bins = int(1+np.log2(pulse.shape[1]))#Officiel de Sturges
count, level = np.histogram(pulse.max(),bins=bins)
level = np.array(list(map(lambda x:(level[x]+level[x+1])/2,range(level.size-1))))
plt.plot(level,count,marker='o',linestyle='None')
plt.show()

Quelque chose comme une distribution normale est apparu.

Montage avec lmfit

Comme vous pouvez le voir depuis import ~, lmfit utilise des objets de classe pour l'ajustement. La fonction d'ajustement est une distribution gaussienne et est la suivante. Le centre du pic est donné par $ \ mu $, et la variation de distribution est donnée par $ \ sigma . amp est le terme d'amplitude. $f_{\left(x\right)}=\frac{amp}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{\left(x-\mu \right)^2}{2\sigma^2}\right)$$ Tout d'abord, donnez une fonction qui correspond à la classe Model. Ensuite, créez une classe Parameters et donnez-lui des valeurs initiales. La classe Model a une méthode .make_params (), qui crée automatiquement un objet Parameters après le deuxième argument de la fonction donnée.

.ipynb

from lmfit import Model, Parameters, Parameter
def gaussian(x,amp,mu,sigma):
    out = amp / (sigma * np.sqrt(2 * np.pi)) * np.exp(-(x - mu)**2 /(2 * sigma**2))
    return out
model = Model(gaussian)#Attribuer une fonction d'ajustement, créer un objet de modèle
pars = model.make_params()#Génération d'objets de paramètres
pars

L'objet Parameters créé est un ensemble d'objets Parameter. Chaque paramètre est accessible comme un dictionnaire. De là avec la méthode .set ・ Valeur de la valeur initiale ・ Limites supérieure et inférieure des paramètres, max, min: -np.inf ~ np.inf -Si pour changer la valeur avant et après ajustement (contrainte) varie = Vrai ou Faux Vous pouvez jouer avec chacun d'eux.

.ipynb

pars['mu'].set(value=70)
pars['amp'].set(value=1)
pars['sigma'].set(value=1)
pars

C'est très approximatif, mais j'ai donné la valeur initiale. Les paramètres peuvent être définis en détail et ils sont affichés dans un tableau de type DataFrame. C'est cool ...

Effectuer le montage

.ipynb

result = model.fit(data=count,params=pars,x=level)
result

Le résultat est affiché comme ceci. Il semble qu'un nouvel objet ModelResult soit créé lorsque .fit est exécuté. L'objet généré avait une méthode .plot_fit (). Visualisez avec cette fonction.

Vous pouvez récupérer les paramètres optimisés dans un format de dictionnaire à l'aide de la méthode .best_values.

.ipynb

opt = result.best_values
opt['mu']/opt['sigma']

La netteté a été évaluée par $ \ frac {\ mu} {\ sigma} $. C'était 23,67.

J'ai essayé de changer la taille du pas sur l'axe horizontal

・ Officiel de Scott $h = \frac{3.49\sigma}{ n^\frac{1}{3} }$ J'ai essayé de définir la taille des pas avec. Comme c'est un gros problème, j'utiliserai celui que j'ai trouvé dans le raccord pour $ \ sigma $.

.ipynb

h = 3.49 * opt['sigma'] / (pulse.shape[1] **(1/3))#Officiel de Scott
r = pulse.max().max()-pulse.max().min()
bins = int(r/h)

count, level = np.histogram(pulse.max(),bins=bins)
level = np.array(list(map(lambda x:(level[x]+level[x+1])/2,range(level.size-1))))
plt.plot(level,count,marker='o',linestyle='None')
plt.show()

La taille des pas est devenue plus étroite et la valeur maximale de la fraction est devenue plus petite, mais la forme de la distribution semble être correcte. Ajuster en échangeant uniquement les valeurs de données et x à substituer. La netteté est de ... 24,44. Améliorée de 0,77.

J'ai essayé de pondérer l'erreur

C'est le cas d'ici.

Il voulait renforcer la plage de $ \ mu \ pm \ sigma $ et ajuster les poids inférieurs à 63 à 0. Il semble que la pondération puisse être effectuée en spécifiant l'option poids dans .fit.

.ipynb

s = (level>opt['mu']-opt['sigma']) * (level<opt['mu']+opt['sigma'])
weight = (np.ones(level.size) + 4*s)*(level>63)
result3 = model.fit(data=count,params=pars,x=level,weights=weight)

↓ Poids et résultat Cela n'a pas l'air très différent, mais si vous regardez de près, vous avez l'impression d'être un peu abaissé. La netteté est de ... 25,50. Seulement 0,06 est devenu plus net.

Les points de post-ajustement peuvent être appelés en tant que type de tableau à l'aide de la méthode .best_fit. C'était pratique lors de la création de csv. Résolvez le problème en toute sécurité ...

Résumé

Cette fois, les paramètres ont été générés automatiquement avec Model.make_params (), mais il existe également une méthode pour générer directement une classe vide avec Parameters () et y ajouter des variables avec la méthode .add.

J'étais confus au début car j'ai utilisé plusieurs objets pour l'ajustement. Il était facile de ne pas savoir où et quelle méthode correspond. Cependant, au fur et à mesure que vous vous y habituerez, vous apprécierez le (s) paramètre (s).

Le fait qu'il y ait peu de sites expliquant lmfit en japonais est aussi la raison pour laquelle j'ai commencé à écrire cet article. Y a-t-il une meilleure bibliothèque ...? Je vous serais reconnaissant si vous pouviez m'apprendre. À l'avenir, les juniors ne seront pas en difficulté. C'était la première fois que j'écrivais du code en lisant la référence officielle (rires). Il n'y a pas d'autre choix que de le faire en raison des résultats de la recherche.

À l'avenir, j'aimerais en apprendre un peu plus sur la technologie qui semble être connectée à ma vie, et je veux également pouvoir utiliser la technologie que j'ai dans divers domaines. Je veux être créatif.

C'est tout pour les progrès liés à python. Je n'ai pas pu m'empêcher de le signaler au professeur (personne ne fait du python), alors j'aimerais profiter de cette occasion pour le signaler. Je suis profondément reconnaissant au professeur d'avoir attendu patiemment l'analyse de mes données.