[PYTHON] Algorithme étonnant "BEADS" qui peut bien séparer l'arrière-plan gênant des données de mesure sans définition de fonction

Récemment, les possibilités d'entrer en contact avec différents types de données de mesure se sont multipliées et le sentiment que «je veux effectuer une évaluation de fond dans un lot sans dépendre de l'échantillon ou de la méthode de mesure…» s'est accru. Après avoir recherché diverses choses, j'ai trouvé que "vous pouvez faire un ajustement en arrière-plan sans spécifier une plage qui évite les pics avec juste un peu d'ingéniosité" Paper viewdoc / download? doi = 10.1.1.64.4698 & rep = rep1 & type = pdf) a été trouvé. Puisqu'il y avait des données de mesure continues pour lesquelles il était difficile de spécifier la plage simplement parce que le pic se déplaçait, je pensais que c'était le mieux que je n'avais pas besoin de spécifier la plage, alors je l'ai immédiatement adoptée et j'ai réalisé l'effet, mais "Peut-être une technologie d'évaluation d'arrière-plan Je creuse un peu plus en pensant "Est-ce que ça progresse?" BEADS J'ai trouvé. L'arrière-plan peut être évalué à grande vitesse et sans modèle (bien qu'il existe des hyper paramètres). Quand je l'ai essayé, j'ai pu le faire, et j'ai pensé: "Cela devrait être plus largement connu ...", alors je vais le présenter.

Contexte

(Si vous voulez voir un exemple rapidement, vous pouvez le lire depuis la préparation)

Suppression d'arrière-plan que je savais

Les données qui sont traitées quotidiennement dans toutes les expériences et mesures telles que l'astronomie, la physique, la chimie et la biologie sont presque toujours appelées l'arrière-plan (base de référence dans certains domaines), et si possible, il y a un signal que je voulais ne pas mesurer. Il contient. Il existe différentes formes dans lesquelles l'arrière-plan est mélangé avec les données de mesure. Par exemple, une forme très simple qui soulève toutes les données de mesure avec une valeur constante ou une ligne droite qui monte / descend vers la droite, ou plusieurs collines douces se chevauchent. Il y a quelque chose qui ressemble. バックグラウンドの例

La première étape de l'analyse des données expérimentales est l'évaluation de ce contexte, mais la méthode d'évaluation réellement utilisée dépend de la méthode expérimentale, de la forme de l'échantillon et des éléments constitutifs, des conditions de mesure, etc. Le plus chanceux est lorsque la cause de l'arrière-plan est connue et peut être calculée théoriquement, ou lorsque les données de fond uniquement peuvent être mesurées. La prochaine chance est que le logiciel fourni par le fabricant de l'équipement de mesure comprend une fonction de suppression automatique de l'arrière-plan, et les données supprimées en arrière-plan apparaissent d'une simple pression sur un bouton. Je pense que la plupart des personnes malchanceuses utilisent probablement la méthode d'ajustement le moins carré en utilisant une fonction (comme un polynôme) qui semble être capable d'exprimer la forme de l'arrière-plan. [^ 1]. La procédure requise à ce moment sera la suivante.

  1. Afficher les données de mesure
  2. Définissez une fonction qui peut exprimer l'arrière-plan (par exemple, combien d'ordres d'un polymorphe doit être utilisé)
  3. Sélectionnez le "point de mesure de fond uniquement qui ne semble pas contenir le vrai signal" à installer.
  4. Trouvez une bonne valeur initiale par essais et erreurs
  5. Fit en utilisant un logiciel ou un package prêt à l'emploi

En fait, 3 "Sélection des points de mesure à utiliser pour l'ajustement" se trouve dans Article introduit au début Il peut être ignoré en utilisant la méthode de rendre la fonction de coût positive, négative et asymétrique, mais cette histoire n'est pas mentionnée ici [^ 2]. Et 4 est toujours un problème avec ou sans 3. Cela peut être gênant, en particulier lorsque vous souhaitez traiter par lots de grandes quantités de données.

Je m'inquiète également de la définition de la fonction. Par exemple, dans l'exemple de "deux" collines "+ constantes" en bas à droite de la figure ci-dessus, les deux grandes fonctions et constantes gaussiennes représentées par les lignes discontinues sont ajoutées au signal vrai. Cependant, peu de gens peuvent juger que «cet arrière-plan peut être exprimé en bloquant constamment deux grandes fonctions gaussiennes» simplement en regardant la ligne bleue. Il peut être possible de l'exprimer raisonnablement bien avec un polynôme qui est une option pour le moment dans de tels cas, mais il semble qu'il sera difficile de déterminer l'ordre et de trouver la valeur initiale, et après tout le traitement par lots comme la mesure continue reste mal à l'aise.

Changements dans la technologie d'évaluation de base depuis 2000

Dans les domaines de la physique physique et des sciences des matériaux que j'ai expérimentés, la correction d'arrière-plan n'est pas considérée comme un volet spécial, et probablement des personnes dont moi-même qui ont adopté la méthode d'ajustement par fonction ou de sélection et de complémentarité des points. Je pense qu'il y en a beaucoup. Cependant, il semble que le mouvement de recherche de meilleures techniques d'évaluation de fond soit actif depuis environ 2000 dans la communauté qui utilise des mesures avec des formes de fond complexes telles que la spectroscopie Raman, l'IF-IR, la chromatographie, l'analyse de masse et la RMN. [^ 3]. Parmi eux, l 'article a annoncé en 2010 l'algorithme appelé airPLS (adaptative itérative repondérée des moindres carrés pénalisés), qui est devenu particulièrement célèbre et est depuis toujours apparu dans les benchmarks. blob / master / airPLS_manuscript.pdf) catégorise les principales méthodes annoncées jusqu'à présent. Selon la méthode générale

Il est divisé en deux, et chacun a ses propres bons et mauvais. Cet article n'a pas l'intention de faire référence aux techniques d'évaluation de base en général, nous ne discuterons donc pas de ces techniques plus loin ici. Si vous êtes intéressé, veuillez suivre le devis airPLS pour le découvrir. Certaines données de mesure peuvent être plus performantes que BEADS.

BEADS BEADS (__B__aseline __E__stimation __A__nd __D__enoising using __S__parsity) a été annoncé en 2014 dans le magazine Chemometrics and Intelligent Laboratory Systems (auteur Duval). Il existe également un [preprint] publié par M. (http://www.laurent-duval.eu/Articles/Ning_X_2014_j-chemometr-intell-lab-syst_chromatogram_bedusbeads-preprint.pdf). Les auteurs sont issus de la communauté du traitement du signal, mais il s'agit d'un algorithme proposé comme solution au problème de traitement de base des chromatographes.

Je ne suis pas très doué pour lire des formules mathématiques, mais le but de cet algorithme est

Je comprends que. Si vous voulez comprendre sur la base de formules mathématiques, veuillez consulter l'article original. Désormais, nous montrerons sa puissance par des exemples concrets.

Préparation

Le BEADS original est écrit en MATLAB, mais [traduction Python] J'ai trouvé (https://github.com/hsiaocy/Beads), et sur cette base, j'ai créé une version qui a amélioré la vitesse d'exécution au même niveau que la version MATLAB. https://github.com/skotaro/pybeads/ Puisqu'il est également publié sur PyPi, il peut être installé en tant que package à l'aide de pip.

pip install pybeads

Cet article utilise les packages suivants, y compris les pybeads.

import numpy as np
import matplotlib.pyplot as plt
import pybeads

Il existe également un notebook Jupyter qui fait la même chose que cet article. https://github.com/skotaro/pybeads/blob/master/sample_data/pybeads_demo.ipynb

Exemple de données

Les exemples de données sont https://github.com/skotaro/pybeads/blob/master/sample_data/chromatograms_and_noise.csv. Il s'agit des données de mesure réelles de la chromatographie attachées au Code MATLAB. / chromatograms.mat et le bruit blanc artificiel data / noise.mat` sont convertis du format MATLAB en fichier CSV et combinés en un seul fichier (l'autorisation de redistribution a été obtenue). Il y a 9 colonnes au total, les 8 premières colonnes sont 8 données réelles avec différents rapports S / N, et les 9 dernières colonnes sont du bruit blanc qui semble avoir été créé pour la démonstration de l'élimination du bruit. Ici, j'utiliserai les 4èmes données avec un fond particulièrement grand plus le bruit [^ 4].

sample_data = np.genfromtxt('sample_data/chromatograms_and_noise.csv', delimiter=',')
y = sample_data[:, 3] + sample_data[:, 8]
print(y.shape)

fig, axes = plt.subplots(2, 1, figsize=(7, 6))
axes[0].plot(y)
axes[1].plot(y)
axes[1].set_ylim(-10, 200)

サンプルデータの様子

Je ne suis pas familier avec la chromatographie, donc je ne sais pas si un tel arrière-plan est courant, mais quand il s'agit d'ajuster ce fond avec un polypole, il semble un peu difficile de trouver la valeur initiale. Appliquons BEADS à ces données et estimons le fond.

fc = 0.006 #Fréquence de coupure utilisée pour créer un filtre passe-haut
d = 1 #Paramètres lors de la création d'un filtre passe-haut. 1 va bien. Consultez l'article pour plus de détails.
r = 6 #Asymétrie de la pénalité d'asymétrie

#Paramètres normalisés pour les données de mesure et leur différenciation
amp = 0.8
lam0 = 0.5 * amp
lam1 = 5 * amp
lam2 = 4 * amp
#Nombre de boucles dans l'algorithme MM
Nit = 15
#Fonction de pénalité
pen = 'L1_v2'

signal_est, bg_est, cost = pybeads.beads(y, d, fc, r, Nit, lam0, lam1, lam2, pen, conv=None)

fig, axes = plt.subplots(3, 1, figsize=(12, 7), sharex=True)
fig.subplots_adjust(hspace=0)
fig.patch.set_color('white')
axes[0].plot(y, c='k', label='original data')
axes[0].plot(bg_est, c='r', label='BG estimated by BEADS')
axes[0].legend()
axes[0].set_ylim(-20, 350)
axes[0].set_xlim(0, 4000)

axes[1].plot(signal_est, label='signal estimated by BEADS')
axes[1].legend()
axes[1].set_ylim(-20, 350)
axes[2].plot(y-signal_est-bg_est, label='noise estimated by BEADS')
axes[2].set_ylim(-35, 35)
axes[2].legend()

BEADSの結果

La ligne noire dans le graphique supérieur représente les données d'exemple. La ligne rouge est le résultat de l'estimation. Êtes-vous surpris? J'ai été honnêtement surpris que "Eh ... c'est incroyable ...". Le temps d'exécution est d'environ 400ms avec 4000 points. Dans l'article utilisant le code MATLAB, 1000 points correspondent à 120 ms, il semble donc que la vitesse soit la même même dans la version Python. De plus, même si je n'y ai pas encore parlé, en fait, comme le montre le graphique en bas, la séparation du bruit est également effectuée en même temps.

Puisque BEADS est basé sur l'algorithme de Majorisation-Minimisation qui fait tourner une boucle et s'approche de la solution, le traçage de la façon dont le coût calculé à partir des données de mesure et de la différenciation (une mesure d'estimation bonne ou mauvaise) est réduit pour chaque boucle est tracé à la solution avec les paramètres définis Vous pouvez voir comment est la convergence. Il semble bon de juger que les données cette fois ont suffisamment convergé dans les 15 temps spécifiés au moment de l'exécution.

コストの遷移

Différence par rapport aux données réelles

En fait, cet exemple de données a une propriété qui convient à BEADS selon laquelle "les valeurs mesurées tombent progressivement à 0 aux deux extrémités des données". Essayez d'ajouter des constantes afin que les deux extrémités soient différentes de zéro et essayez d'estimer en utilisant les mêmes paramètres.

bg_const = 50
signal_est, bg_est, cost = pybeads.beads(y+bg_const, d, fc, r, Nit, lam0, lam1, lam2, pen, conv=None)

両端が0でなくなるとダメ

Je ne peux pas bien estimer aux deux extrémités. Vous pouvez également voir que l'estimation est de 0 aux deux extrémités. Dans les données de mesure réelles, la valeur mesurée ne devient pas toujours égale à 0 aux deux extrémités. C'est un problème. Que faire maintenant?

S'il ne tombe pas (à 0), laissez-le tomber. Valeur mesurée

Je l'ai essayé avec les données à portée de main et cela n'a pas fonctionné, et j'ai été surpris de constater que "c'est un algorithme tellement merveilleux, mais ce doit être des données de mesure qui tombent en douceur à 0 aux deux extrémités ...", mais j'ai remarqué l'astuce de données d'échantillonnage ci-dessus et en douceur N'est-il pas possible d'étendre les données de mesure en créant une pièce qui devient 0? Il a clignoté [^ 5].

Tout d'abord, essayons de rendre l'arrière-plan plus suspect que d'élever le fond avec une constante.

bg = 5e-5*(np.linspace(0, 3999, num=4000)-2000)**2
y_difficult = y + bg
plt.plot(y_difficult)
plt.ylim(0, 350)

もっとやらしいバックグラウンド

C'est devenu très simple. Ajoutez une extension à cela. La fonction sigmoïde est utilisée pour la fonction qui tombe du lécher. xscale_l etc. sont des paramètres qui" étendent "la partie de la fonction sigmoïde qui tombe à 0, et plus le nombre est grand, plus l'extension" étend ". En d'autres termes, de nombreux pseudo points de mesure sont ajoutés à la partie d'extension avant qu'elle ne tombe à 0. Pour le moment, définissons le paramètre d'extension sur 30 pour la gauche et la droite.

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

xscale_l, xscale_r = 30, 30
dx = 1
y_difficult_l = y_difficult[0]*sigmoid(1/xscale_l*np.arange(-5*xscale_l, 5*xscale_l, dx))
y_difficult_r = y_difficult[-1]*sigmoid(-1/xscale_r*np.arange(-5*xscale_r, 5*xscale_r, dx))
y_difficult_ext = np.hstack([y_difficult_l, y_difficult, y_difficult_r])
len_l, len_o, len_r = len(y_difficult_l), len(y_difficult), len(y_difficult_r)
plt.plot(range(len_l, len_l+len_o), y_difficult)
plt.plot(y_difficult_l, 'C1')
plt.plot(range(len_l+len_o, len_l+len_o+len_r), y_difficult_r, 'C1')

両端を延長させた

La ligne orange est l'extension. Cela semble bon d'une certaine manière, alors je vais l'estimer.

signal_est, bg_est, cost = pybeads.beads(y_difficult_ext, d, fc, r, Nit, lam0, lam1, lam2, pen, conv=None)

惜しい

Cela ne semble pas être retardé, mais cela semble être une bonne approche. Définissons le paramètre d'extension sur 100.

xscale_l, xscale_r = 100, 100

もっと間延びさせてみる

signal_est, bg_est, cost = pybeads.beads(y_difficult_ext, d, fc, r, Nit, lam0, lam1, lam2, pen, conv=None)

うまくいった！

Réussi! Cela traitera probablement toutes les données de mesure (si le pic est quelque peu clairsemé).

À propos de l'ajustement de para élevé

En bref, BEADS est sans modèle et peut estimer l'arrière-plan sans définir de fonction, mais il a tellement d'hyperparamètres.

fc = 0.006 #Fréquence de coupure utilisée pour créer un filtre passe-haut
d = 1 #Paramètres lors de la création d'un filtre passe-haut. 1 va bien. Consultez l'article pour plus de détails.
r = 6 #Asymétrie de la pénalité d'asymétrie

#Paramètres normalisés pour les données de mesure et leur différenciation
amp = 0.8
lam0 = 0.5 * amp
lam1 = 5 * amp
lam2 = 4 * amp
#Nombre de boucles dans l'algorithme MM
Nit = 15
#Fonction de pénalité
pen = 'L1_v2'

Je ne suis pas encore habitué à la sensation de la taille des nombres, mais si je joue un peu avec, le plus important est fc, qui décide quelle fréquence ou moins le signal doit être en arrière-plan, et la suivante est efficace. Paramètre de normalisation. En fait, en fonction des données, même une petite modification du paramètre de normalisation (par exemple, 0,001 à 0,0011) peut faire une grande différence dans le résultat. Dans ce cas, essayez conv = 3 ou conv = 5. Le résultat différentiel est lissé par la moyenne mobile et le résultat est stable [^ 6].

En ce qui concerne d, r, Nit, pen, je pense que cela dépend des données, mais je n'ai pas eu besoin de le changer avec les données que je gère.

Résumé

--BEADS permet une estimation sans modèle et la séparation du bruit de fond et du bruit des données de mesure sans définir de fonction.

Matériel de référence

[^ 1]: Il existe également une méthode primitive mais fiable de complément de spline ou de complément linéaire des points sélectionnés. Ce n'est pas une mauvaise méthode lorsque la quantité de données que vous souhaitez traiter n'est pas si grande ou lorsque le flou causé par l'opération manuelle ne vous dérange pas. J'évite cette méthode autant que possible car il y a beaucoup de données de mesure que je souhaite traiter et j'ai peur du flou lorsque je le fais manuellement, et surtout cela me rend triste. Il existe également un moyen d'utiliser la conversion en ondelettes, mais avant de l'examiner en détail, j'étais désintéressé car il est devenu "Je suis d'accord avec BEADS".

[^ 2]: Si vous connaissez une fonction appropriée pour l'évaluation d'arrière-plan, ou si vous devez évaluer par fonction, la méthode de la fonction de coût asymétrique est assez efficace et vaut la peine d'être essayée. Il existe également une boîte à outils matlab de l'auteur de l'article, mais vous pouvez facilement l'écrire en Python en utilisant scipy.optimize.minimize. ..

[^ 3]: Étant donné que ces méthodes de mesure sont largement utilisées dans des domaines avec de nombreux chercheurs de l'industrie et du monde universitaire tels que la chimie, la biologie, la médecine et l'alimentation, il y avait tout simplement beaucoup de demande et l'impact de l'amélioration était grand, ce qui est l'une des raisons pour lesquelles l'amélioration a progressé. Je me demande si.

[^ 4]: En regardant la valeur de sample_data, il y a une valeur négative avant même d'ajouter du bruit blanc. Peut-être que les données réelles ont déjà été traitées pour tracer une ligne reliant les deux extrémités des données de mesure brutes.

[^ 5]: Cela ne falsifie aucune donnée de mesure. Il n'y a aucun problème si vous l'étendez uniquement pour estimer l'arrière-plan de la partie de mesure et supprimez la partie d'extension après avoir soustrait l'arrière-plan.

[^ 6]: C'est une fonctionnalité que j'ai ajoutée, donc elle n'est pas incluse dans la version originale de MATLAB (pour l'instant).

Recommended Posts

Algorithme étonnant "BEADS" qui peut bien séparer l'arrière-plan gênant des données de mesure sans définition de fonction