[PYTHON] Erstaunlicher Algorithmus "BEADS", der den störenden Hintergrund von Messdaten ohne Funktionsdefinition gut trennen kann

In letzter Zeit haben die Möglichkeiten, mit verschiedenen Arten von Messdaten in Kontakt zu kommen, zugenommen, und das Gefühl, dass "ich eine Hintergrundbewertung in einer Charge durchführen möchte, ohne von der Probe oder der Messmethode abhängig zu sein ...", hat zugenommen. Nachdem ich nach verschiedenen Dingen gesucht hatte, stellte ich fest, dass "Sie eine Hintergrundanpassung durchführen können, ohne einen Bereich anzugeben, der Spitzen mit nur wenig Einfallsreichtum vermeidet" Paper viewdoc / download? doi = 10.1.1.64.4698 & rep = rep1 & type = pdf) wurde gefunden. Da es kontinuierliche Messdaten gab, bei denen es schwierig war, den Bereich zu spezifizieren, nur weil sich der Peak bewegte, dachte ich, dass es das Beste war, dass ich den Bereich nicht spezifizieren musste, also übernahm ich ihn sofort und erkannte den Effekt, aber "Vielleicht Hintergrundbewertungstechnologie Ich grabe etwas tiefer und denke: "Geht es voran?" BEADS Ich fand. Der Hintergrund kann mit hoher Geschwindigkeit und ohne Modell ausgewertet werden (obwohl es Hyperparameter gibt). Als ich es versuchte, war ich in der Lage, es zu tun, und ich dachte, "Das sollte bekannter sein ...", also werde ich es vorstellen.

Hintergrund

(Wenn Sie ein Beispiel schnell sehen möchten, können Sie es aus der Vorbereitung lesen.)

Hintergrundentfernung wusste ich

Die Daten, die täglich in allen Experimenten und Messungen wie Astronomie, Physik, Chemie und Biologie verarbeitet werden, werden fast immer als Hintergrund bezeichnet (Basislinie in einigen Bereichen), und wenn möglich, gibt es ein Signal, dass ich nicht gemessen werden wollte. Es beinhaltet. Es gibt verschiedene Formen, in denen der Hintergrund mit den Messdaten gemischt wird. Zum Beispiel eine sehr einfache Form, bei der die gesamten Messdaten mit einem konstanten Wert oder einer geraden Linie, die nach rechts steigt / fällt, oder mehrere sanfte Hügel überlappen. Es gibt etwas, das aussieht. バックグラウンドの例

Der erste Schritt bei der Analyse experimenteller Daten ist die Bewertung dieses Hintergrunds. Die tatsächlich verwendete Bewertungsmethode hängt jedoch von der experimentellen Methode, der Probenform und den Bestandteilen, den Messbedingungen usw. ab. Am glücklichsten ist es, wenn die Ursache des Hintergrunds bekannt ist und theoretisch berechnet werden kann oder wenn die Nur-Hintergrund-Daten gemessen werden können. Das nächste Glück ist, dass die vom Hersteller der Messgeräte bereitgestellte Software eine automatische Hintergrundentfernungsfunktion enthält und die vom Hintergrund entfernten Daten per Knopfdruck angezeigt werden. Ich denke, die meisten unglücklichen Menschen verwenden wahrscheinlich die Methode der kleinsten Quadrate mithilfe einer Funktion (z. B. eines Polynoms), die die Form des Hintergrunds ausdrücken kann. [^ 1]. Das zu diesem Zeitpunkt erforderliche Verfahren ist wie folgt.

  1. Zeigen Sie die Messdaten an
  2. Definieren Sie eine Funktion, die den Hintergrund ausdrücken kann (z. B. wie viele Ordnungen eines Polymorphs verwendet werden sollen).
  3. Wählen Sie den "Nur-Hintergrund-Messpunkt, der nicht das wahre Signal zu enthalten scheint", der angepasst werden soll.
  4. Finden Sie durch Ausprobieren einen guten Anfangswert
  5. Mit einer vorgefertigten Software oder einem Paket montieren

Tatsächlich steht 3 "Auswahl der für die Anpassung zu verwendenden Messpunkte" in Artikel zu Beginn eingeführt Es kann übersprungen werden, indem die Kostenfunktion positiv und negativ und asymmetrisch gemacht wird, aber diese Geschichte wird hier nicht erwähnt [^ 2]. Und 4 ist immer noch ein Ärger mit oder ohne 3. Dies kann unpraktisch sein, insbesondere wenn Sie große Datenmengen stapelweise verarbeiten möchten.

Ich mache mir auch Sorgen um die Definition der Funktion. Zum Beispiel werden im Beispiel von "zwei" Hügeln "+ Konstanten" unten rechts in der obigen Abbildung die zwei breiten Gaußschen Funktionen und Konstanten, die durch die gestrichelten Linien dargestellt sind, zum wahren Signal addiert. Nur wenige Menschen können jedoch beurteilen, dass "dieser Hintergrund ausgedrückt werden kann, indem zwei breite Gaußsche Funktionen ständig verstopft werden", wenn man nur die blaue Linie betrachtet. Es mag möglich sein, es mit einem Polynom, das in solchen Fällen vorerst eine Option ist, einigermaßen gut auszudrücken, aber es scheint schwierig zu sein, die Reihenfolge zu bestimmen und den Anfangswert zu finden, und schließlich bleibt die Stapelverarbeitung wie die kontinuierliche Messung unruhig.

Änderungen in der Hintergrundbewertungstechnologie ab etwa 2000

In den Bereichen der physikalischen Physik und Materialwissenschaften, die ich erlebt habe, wird die Hintergrundkorrektur nicht als spezieller Bereich angesehen, und wahrscheinlich auch für Personen, die die Methode der Anpassung nach Funktion oder Auswahl und Ergänzung von Punkten übernommen haben. Ich denke es gibt viele. Es scheint jedoch, dass die Bewegung zur Suche nach besseren Hintergrundbewertungstechniken seit etwa 2000 in der Gemeinde aktiv ist, die Messungen mit komplizierten Hintergrundformen wie Raman-Spektroskopie, IF-IR, Chromatographie, Massenanalyse und NMR verwendet. [^ 3]. Darunter der 2010 veröffentlichte Artikel des Algorithmus airPLS (adaptive iterativ neu gewichtete bestrafte kleinste Quadrate), der besonders berühmt wurde und seitdem immer in den Benchmarks auftauchte. blob / master / airPLS_manuscript.pdf) kategorisiert die bisher angekündigten Hauptmethoden. Nach der allgemeinen Methode

--Polygonanpassung

Es ist zweigeteilt und jedes hat sein eigenes Gut und Böse. Dieser Artikel beabsichtigt nicht, sich allgemein auf Hintergrundbewertungstechniken zu beziehen, daher werden wir diese Techniken hier nicht weiter diskutieren. Wenn Sie interessiert sind, folgen Sie bitte dem airPLS-Angebot, um dies herauszufinden. Einige Messdaten sind möglicherweise besser als BEADS.

BEADS BEADS (__B__aseline __E__stimation __A__nd __D__enoising using __S__parsity) wurde 2014 in der Zeitschrift Chemometrics and Intelligent Laboratory Systems (Autor Duval) angekündigt. Es gibt auch einen [Preprint], der von Mr. veröffentlicht wurde (http://www.laurent-duval.eu/Articles/Ning_X_2014_j-chemometr-intell-lab-syst_chromatogram_bedusbeads-preprint.pdf). Die Autoren stammen aus der Signalverarbeitungsgemeinschaft, dies ist jedoch ein Algorithmus, der als Lösung für das Grundverarbeitungsproblem von Chromatographen vorgeschlagen wird.

Ich bin nicht sehr gut darin, mathematische Formeln zu lesen, aber der Sinn dieses Algorithmus ist

Ich verstehe das. Wenn Sie anhand mathematischer Formeln verstehen möchten, lesen Sie bitte das Originalpapier. Von nun an werden wir seine Kraft anhand konkreter Beispiele zeigen.

Vorbereitung

Das Original BEADS ist in MATLAB geschrieben, aber [Python-Übersetzung] Ich fand (https://github.com/hsiaocy/Beads) und erstellte darauf basierend eine Version, die die Ausführungsgeschwindigkeit auf das gleiche Niveau wie die MATLAB-Version verbesserte. https://github.com/skotaro/pybeads/ Da es auch auf PyPi veröffentlicht wird, kann es mit pip als Paket installiert werden.

pip install pybeads

In diesem Artikel werden die folgenden Pakete verwendet, einschließlich Pybeads.

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

Es gibt auch ein Jupyter-Notizbuch, das das Gleiche wie dieser Artikel tut. https://github.com/skotaro/pybeads/blob/master/sample_data/pybeads_demo.ipynb

Beispieldaten

Die Beispieldaten sind https://github.com/skotaro/pybeads/blob/master/sample_data/chromatograms_and_noise.csv. Dies sind die tatsächlichen Messdaten der Chromatographie, die an [MATLAB-Code] angehängt sind (https://jp.mathworks.com/matlabcentral/fileexchange/49974-beads-baseline-estimation-and-denoising-with-sparsity). / chromatograms.matund künstliches weißes Rauschendata / Noise.mat` werden vom MATLAB-Format in eine CSV-Datei konvertiert und zu einer Datei zusammengefasst (Weiterverteilungserlaubnis wurde eingeholt). Insgesamt gibt es 9 Spalten, die ersten 8 Spalten sind 8 tatsächliche Daten mit unterschiedlichen S / N-Verhältnissen, und die letzten 9 Spalten sind weißes Rauschen, das anscheinend zur Demonstration der Rauschentfernung erstellt wurde. Hier werde ich die 4. Daten mit einem besonders großen Hintergrund plus Rauschen verwenden [^ 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)

サンプルデータの様子

Ich bin mit Chromatographie nicht vertraut, daher weiß ich nicht, ob ein solcher Hintergrund üblich ist, aber wenn es darum geht, diesen Hintergrund mit einem Polypoly zu versehen, scheint es ein wenig schwierig zu sein, den Anfangswert zu finden. Wenden wir BEADS auf diese Daten an und schätzen den Hintergrund.

fc = 0.006 #Grenzfrequenz zur Erstellung des Hochpassfilters
d = 1 #Parameter beim Erstellen eines Hochpassfilters. 1 ist in Ordnung. Einzelheiten finden Sie im Papier.
r = 6 #Asymmetrie der Asymmetriestrafe

#Normalisierte Parameter für Messdaten und deren Differenzierung
amp = 0.8
lam0 = 0.5 * amp
lam1 = 5 * amp
lam2 = 4 * amp
#Anzahl der Schleifen im MM-Algorithmus
Nit = 15
#Straffunktion
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の結果

Die schwarze Linie im oberen Diagramm zeigt die Beispieldaten. Die rote Linie ist das Schätzergebnis. Bist du überrascht? Ich war ehrlich überrascht, dass "Eh ... es ist erstaunlich ...". Die Ausführungszeit beträgt ca. 400ms mit 4000 Punkten. In der Arbeit mit MATLAB-Code sind 1000 Punkte 120 ms, daher scheint die Geschwindigkeit auch in der Python-Version gleich zu sein. Auch wenn ich es bisher noch nicht angesprochen habe, wird, wie in der Grafik unten gezeigt, gleichzeitig auch die Rauschtrennung durchgeführt.

Da BEADS auf dem Majorization-Minimization-Algorithmus basiert, der eine Schleife dreht und sich der Lösung nähert, wird die Lösung aus den Messdaten und der Differenzierung (ein Maß für die gute oder schlechte Schätzung) für jede Schleife mit den eingestellten Parametern in die Lösung aufgetragen. Sie können sehen, wie die Konvergenz ist. Es scheint gut zu beurteilen, dass die Daten dieses Mal in den zum Zeitpunkt der Ausführung angegebenen 15 Zeiten ausreichend konvergiert haben.

コストの遷移

Unterschied zu tatsächlichen Daten

Tatsächlich haben diese Beispieldaten eine für BEADS geeignete Eigenschaft, dass "die gemessenen Werte an beiden Enden der Daten allmählich auf 0 fallen". Versuchen Sie, Konstanten so hinzuzufügen, dass beide Enden ungleich Null sind, und versuchen Sie, mit denselben Parametern zu schätzen.

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

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

Ich kann an beiden Enden nicht gut einschätzen. Sie können auch sehen, dass die Schätzung an beiden Enden 0 ist. In tatsächlichen Messdaten wird der gemessene Wert an beiden Enden nicht immer gleichmäßig 0. Das ist ein Problem. Was nun?

Wenn es nicht fällt (auf 0), lassen Sie es fallen. Messwert

Ich habe es mit den vorliegenden Daten versucht und es hat nicht funktioniert, und ich war überrascht, dass "es ein so wunderbarer Algorithmus ist, aber es müssen Messdaten sein, die an beiden Enden reibungslos auf 0 fallen ...", aber ich habe den obigen Trick mit den Beispieldaten und reibungslos bemerkt Ist es nicht möglich, die Messdaten zu erweitern, indem ein Teil erstellt wird, das 0 wird? Es blitzte [^ 5].

Lassen Sie uns zunächst versuchen, den Hintergrund verdächtiger zu gestalten, als den Boden mit einer Konstanten anzuheben.

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

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

Es ist sehr einfach geworden. Fügen Sie diesem eine Erweiterung hinzu. Die Sigmoid-Funktion wird für die Funktion verwendet, die vom Leck fällt. "xscale_l" usw. sind Parameter, die den Teil der Sigmoid-Funktion "erweitern", der auf 0 fällt, und je größer die Zahl, desto "erweitern" Sie die Erweiterung. Mit anderen Worten, dem Erweiterungsteil werden viele Pseudomesspunkte hinzugefügt, bevor er auf 0 fällt. Lassen Sie uns vorerst den Erweiterungsparameter für links und rechts auf 30 setzen.

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')

両端を延長させた

Die orange Linie ist die Erweiterung. Es sieht irgendwie gut aus, also werde ich es schätzen.

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

惜しい

Es scheint nicht verzögert zu sein, aber es sieht nach einem schönen Ansatz aus. Setzen wir den Erweiterungsparameter auf 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)

うまくいった！

Erfolgreich! Dies wird wahrscheinlich alle Messdaten verarbeiten (wenn der Peak etwas dünn ist).

Über hohe Para-Anpassung

BEADS ist modellfrei und kann den Hintergrund schätzen, ohne eine Funktion zu definieren, hat aber so viele Hyperparameter.

fc = 0.006 #Grenzfrequenz zur Erstellung des Hochpassfilters
d = 1 #Parameter beim Erstellen eines Hochpassfilters. 1 ist in Ordnung. Einzelheiten finden Sie im Papier.
r = 6 #Asymmetrie der Asymmetriestrafe

#Normalisierte Parameter für Messdaten und deren Differenzierung
amp = 0.8
lam0 = 0.5 * amp
lam1 = 5 * amp
lam2 = 4 * amp
#Anzahl der Schleifen im MM-Algorithmus
Nit = 15
#Straffunktion
pen = 'L1_v2'

Ich bin noch nicht an das Gefühl der Größe von Zahlen gewöhnt, aber wenn ich ein wenig damit spiele, ist das Wichtigste "fc", das entscheidet, welche Frequenz oder weniger das Signal im Hintergrund sein soll, und die nächste ist effektiv. Normalisierungsparameter. Abhängig von den Daten kann sogar eine kleine Änderung des Normalisierungsparameters (z. B. 0,001 bis 0,001) einen großen Unterschied im Ergebnis bewirken. Versuchen Sie in diesem Fall "conv = 3" oder "conv = 5". Das Differentialergebnis wird durch den gleitenden Durchschnitt geglättet und das Ergebnis ist stabil [^ 6].

In Bezug auf "d", "r", "Nit", "Stift" denke ich, dass es von den Daten abhängt, aber ich musste es nicht mit den Daten ändern, die ich verarbeite.

Zusammenfassung

--BEADS ermöglicht die modellfreie Schätzung und Trennung von Hintergrund und Rauschen von Messdaten, ohne eine Funktion zu definieren.

Referenzmaterial

[^ 1]: Es gibt auch eine primitive, aber zuverlässige Methode zur Spline-Komplementierung oder linearen Komplementierung der ausgewählten Punkte. Dies ist keine schlechte Methode, wenn die Datenmenge, die Sie verarbeiten möchten, nicht so groß ist oder wenn Sie die durch manuelle Bedienung verursachten Unschärfen nicht stören. Ich vermeide diese Methode so weit wie möglich, weil ich viele Messdaten verarbeiten möchte und mir Sorgen über Unschärfen bei der manuellen Ausführung mache. Vor allem macht es mich traurig. Es gibt auch eine Möglichkeit, die Wavelet-Konvertierung zu verwenden, aber bevor ich sie im Detail untersuchte, war ich desinteressiert, weil sie zu "Ich bin mit BEADS in Ordnung" wurde.

[^ 2]: Wenn Sie eine geeignete Funktion für die Hintergrundbewertung kennen oder nach Funktionen bewerten müssen, ist die Methode der asymmetrischen Kostenfunktion sehr effektiv und einen Versuch wert. Es gibt auch eine matlab-Toolbox des Autors des Papiers, aber Sie können sie einfach in Python mit "scipy.optimize.minimize" schreiben. ..

[^ 3]: Da diese Messmethoden in Bereichen mit vielen Forschern aus Industrie und Wissenschaft wie Chemie, Biologie, Medizin und Lebensmittel weit verbreitet sind, gab es einfach eine große Nachfrage und die Auswirkungen der Verbesserung waren groß, was ein Grund für den Fortschritt der Verbesserung ist. Ich frage mich, ob.

[^ 4]: Wenn man den Wert von "sample_data" betrachtet, gibt es einen negativen Wert, noch bevor weißes Rauschen hinzugefügt wird. Möglicherweise wurden die tatsächlichen Daten bereits verarbeitet, um eine Linie zwischen den beiden Messdaten zu ziehen, die beide Enden verbindet.

[^ 5]: Hiermit werden keine Messdaten gefälscht. Es ist kein Problem, wenn Sie es nur erweitern, um den Hintergrund des Messteils zu schätzen und den Erweiterungsteil nach dem Subtrahieren des Hintergrunds zu entfernen.

[^ 6]: Dies ist eine Funktion, die ich hinzugefügt habe, daher ist sie (vorerst) nicht in der ursprünglichen MATLAB-Version enthalten.

Recommended Posts

Erstaunlicher Algorithmus "BEADS", der den störenden Hintergrund von Messdaten ohne Funktionsdefinition gut trennen kann