Verwenden Sie das Paketdichtungsverhältnis für das R-Dichteverhältnis von Python

Diesmal handelt es sich bei dem Artikel um eine Kombination der beiden Themen. Zum einen soll das R-Paket von Python aufgerufen und verwendet werden, und zum anderen soll versucht werden, Anomalien mithilfe des Dichteverhältnisses der Dichteverhältnisschätzung R zu erkennen. ist. R ist sehr attraktiv, weil so viele Pakete statistischer Methoden veröffentlicht wurden (mehr als 8000) und einige davon nicht in Python enthalten sind. Wenn Sie jedoch an Python gewöhnt sind, möchten Sie möglicherweise Daten vorverarbeiten und Diagramme in dem bekannten Python zeichnen. (Es geht um mich: heat_smile :) Für diese Leute ist hier, wie man das R-Paket von Python verwendet.

Den vollständigen Text des Python-Quellcodes finden Sie unter [hier] auf GitHub (https://github.com/matsuken92/Qiita_Contents/blob/master/General/Density_ratio_R.ipynb). Verwenden Sie ihn daher entsprechend.

</ i> Umgebung

Die Umgebung, die ich ausprobiert habe, ist wie folgt. Ich versuche es auf Mac und Anaconda. Wenn Sie Anaconda installiert haben, ist keine besondere Vorbereitung erforderlich. Der Inhalt dieses Artikels wird grundsätzlich mit Jupyter Notebook (iPython Notebook) überprüft.

Python 3.5.1 |Anaconda custom (x86_64)| (default, Dec  7 2015, 11:24:55) 
[GCC 4.2.1 (Apple Inc. build 5577)] on darwin

IPython 5.0.0

</ i> R-Paket von Python aus aufrufen

Vorbereitung

Zunächst muss R installiert sein. Starten Sie dann R (oder RStudio) einmal und installieren Sie das Dens-Ratio-Paket wie unten gezeigt. Es wird zur Schätzung des Dichteverhältnisses verwendet. Sobald die Installation abgeschlossen ist, können Sie R beenden.

r


install.packages("densratio")

Greifen Sie über Python auf R zu

Installieren Sie eine Bibliothek namens "rpy2", um von Python aus auf R zuzugreifen.

bash


pip install rpy2

Wir werden mit diesem rpy2 auf R zugreifen. Importieren Sie zunächst verschiedene Dinge. Jetzt können Sie R-Objekte verwenden und zwischen Pandas-Datenrahmen und R-Datenrahmen konvertieren.

Python


from rpy2.robjects.packages import importr
import rpy2.robjects as robjects

# import pandas.rpy.common as com   # [depricated]
from rpy2.robjects import pandas2ri
pandas2ri.activate()

</ i> Probieren Sie das Dichteverhältnis-Schätzpaket für das Dichteverhältnis aus

Ich werde versuchen, das Dichteverhältnis mithilfe des von Boss @hoxo_m erstellten Paketdichteverhältnisses und der Erkennung von Anomalien zu schätzen.

Ich werde versuchen, das R dens-Paket sofort mit Python zu laden, aber ich kann es importieren, indem ich nur eine Zeile schreibe, wie unten gezeigt.

Python


densratio = importr("densratio")

Versuchen Sie es mit eindimensionalen Daten

Die Dichteverhältnisschätzung ist ein Verfahren zum direkten Schätzen des Dichteverhältnisses, anstatt jede Verteilung zu schätzen und dann das Verhältnis zu nehmen, wenn die Verteilung der Eingabevariablen, die als Kovariatenverschiebung bezeichnet wird, zwischen Trainingsdaten und Testdaten unterschiedlich ist. (Eine Kovariate ist eine Eingabevariable)

Trainingsdaten $ \ mathcal {D} = \ {\ boldsymbol {x} ^ {(1)}, \ cdots, \ boldsymbol {x} ^ {(N)} \} $ sind nur normale Daten (wenn abnormal) Selbst wenn es Daten gibt, ist es eine sehr kleine Zahl) und Testdaten $ \ mathcal {D} '= \ {\ boldsymbol {x} ^ {' (1)}, \ cdots, \ boldsymbol {x} ^ {' (N ')} \} $ enthält eine bestimmte Anzahl anomaler Daten, und diese anomalen Daten werden aus dem Dichteverhältnis erkannt.

Wenn die Verteilung der Trainingsdaten $ p_ {tr} (\ boldsymbol {x}) $ und die Verteilung der Testdaten $ p_ {te} (\ boldsymbol {x}) $ ist, wie hoch ist das Dichteverhältnis?

w(\boldsymbol{x}) = {p_{tr}(\boldsymbol{x}) \over p_{te}(\boldsymbol{x})}

Es ist ein Wert ausgedrückt als. Schätzen Sie $ \ hat {w} (\ boldsymbol {x}) $ für dieses $ w (\ boldsymbol {x}) $

\hat{w}(\boldsymbol{x}) = \sum_{l=1}^{N} \alpha_l \varphi_l (\boldsymbol{x})

Und es wird durch eine lineare Kombination der Basisfunktion $ \ {\ varphi_l (\ boldsymbol {x}) \} _ {l = 1} ^ {N} $ dargestellt. Hier wird $ \ varphi_l (\ boldsymbol {x}) $ als Basisfunktion bezeichnet, und diesmal wird der RBF-Kernel für diese Basisfunktion verwendet.

\varphi_l (\boldsymbol{x}) = \exp \left( { \| \boldsymbol{x} - \boldsymbol{x}^{l}\|^2 \over 2\sigma^2} \right)

Wird genutzt. Da $ \ boldsymbol {x} ^ {l} $ die $ l $ -ten Trainingsdaten darstellt, sieht es wie eine Normalverteilung aus, wobei die Normalisierungskonstanten, die auf den Trainingsdaten zentriert sind, weggelassen werden.

Das Folgende ist ein Beispiel unter der Annahme, dass es 5 Trainingsdaten "[0,5, 4,5, 5,0, 7,0, 9,0]" gibt, von denen "α = [0,1, 0,25, 0,35, 0,2, 0,1]" Die gewichteten Kernel werden addiert und der untere Teil zeigt, wie es als geschätztes Dichteverhältnis $ \ hat {w} (\ boldsymbol {x}) $ dargestellt wird.

img012.png

Bei dieser Schätzung des Dichteverhältnisses geht es also darum, das Kernel-Mix-Verhältnis $ \ boldsymbol {\ alpha} $ zu bestimmen. $ \ Sigma $, das die Breite des Kernels angibt, ist ein Hyperparameter, der durch Kreuzvalidierung bestimmt wird.

In den oben genannten Büchern, Seiten und Artikeln erfahren Sie, warum dies funktioniert und wie Sie $ \ boldsymbol {\ alpha} $ entscheiden.

Erzeugung künstlicher Daten

Das übliche Ausreißerproblem besteht darin, den Anomaliegrad einzeln zu berechnen und mit dem Schwellenwert zu vergleichen, unter der Annahme, dass eine neue Daten erhalten werden. Diese Methode ist jedoch dadurch gekennzeichnet, dass es sich bei den Testdaten um mehrere Stichproben handelt. .. Erstens werden diese Daten künstlich unter Verwendung von Zufallszahlen erzeugt.

Da Densratio ein R-Paket ist, verwenden Sie "pandas2ri.py2ri", um Daten von der Pandas-Welt in die R-Welt zu übertragen.

Python


rd.seed(71)
n_data1 = 3000
n_data2 = 500
m = [1, 1.2]
sd = [1/8, 1/4]

x = rd.normal(loc=m[0], scale=sd[0], size=n_data1)
y = rd.normal(loc=m[1], scale=sd[1], size=n_data2)

#Konvertieren Sie Python-Datenrahmen in R-Datenrahmen!
r_x = pandas2ri.py2ri(pd.DataFrame(x))
r_y = pandas2ri.py2ri(pd.DataFrame(y))

Abnormale Daten haben einen etwas höheren Durchschnitt und eine größere Varianz als normale Daten.

img001.png

Führen Sie eine Schätzung des Dichteverhältnisses durch

Berechnen Sie das Dichteverhältnis unter Verwendung des Dichteverhältnisses für diese Daten. In diesem Paket können zwei Arten von Schätzverfahren für das Dichteverhältnis verwendet werden: KLIEP, das das Dichteverhältnis unter Verwendung der Menge an Calback Libra-Informationen schätzt, und uLSIF, das das uneingeschränkte Schätzverfahren für das minimale quadratische Dichteverhältnis verwendet. Dieses Mal verwenden wir uLSIF, das ebenfalls die Standardeinstellung ist und Hochgeschwindigkeitsberechnungen durchführen kann. Fügen Sie einfach das in einen R-Datenrahmen konvertierte r_x, r_y in diese Funktion ein. Es ist praktisch: grinsen: (Dieses Mal wurde $ \ sigma $ ohne Kreuzvalidierung auf 0,2 gesetzt.)

Python


result = densratio.densratio(r_x, r_y, "uLSIF",0.2, "auto")

Während dieses Vorgangs wird eine Warnmeldung angezeigt, die jedoch ein Protokoll des Prozesses zum Bestimmen von Hyperparametern durch Kreuzvalidierung anstelle von Warnung ist. In den rpy2-Spezifikationen wird sie jedoch als Warnung angezeigt.

out


//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: ################## Start uLSIF ##################

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: Searching optimal sigma and lambda...

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.200, lambda = 0.001, score = 11.113

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.200, lambda = 0.003, score = 1.026

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.200, lambda = 0.010, score = -0.654

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.200, lambda = 0.032, score = -0.931

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.200, lambda = 0.100, score = -1.009

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.200, lambda = 0.316, score = -1.029

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: Found optimal sigma = 0.200, lambda = 0.316.

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: Optimizing alpha...

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: End.

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: ################## Finished uLSIF ###############

  warnings.warn(x, RRuntimeWarning)

Extrahieren und führen Sie compute_density_ratio aus, eine Funktion zum Extrahieren des Dichteverhältnisses aus dem Ergebnis auf der Python-Seite. Ermitteln Sie zunächst das geschätzte Dichteverhältnis "w_hat" zur normalen Probe. Wenn das Dichteverhältnis $ \ hat {w} $ ist, ist der Grad der Anomalie minus dem Logarithmus.

a(\boldsymbol{x}^{'(n)}) = -\ln \hat{w}(\boldsymbol{x}^{'(n)})

Da dies durch ausgedrückt wird, wird diese Umwandlung durchgeführt und berechnet.

Darüber hinaus möchte ich eine Abnormalitätsbeurteilung vornehmen, bei der das 99% -Perzentil der Verteilung dieser Trainingsdaten als Abnormalitätsschwelle verwendet wird.

Python


#Berechnen und extrahieren Sie das Dichteverhältnis
compute_density_ratio = get_method(result, "compute_density_ratio")
w_hat = compute_density_ratio(r_x)

#Dichteverhältnis in Anomalie umwandeln
anom_score = -np.log(w_hat)
thresh = np.percentile(anom_score, 99)    #99 des Dichteverhältnisses normaler Proben%Legen Sie das Perzentil als Schwellenwert für die Beurteilung von Anomalien fest

Unten ist eine Darstellung dieser w_hat und anom_score.

img002.png

Berechnen und zeichnen Sie als nächstes das Dichteverhältnis der Testdaten und das Histogramm des Grads der Anomalie.

Python


#Berechnen Sie aus dem Ergebnis_density_Nehmen Sie das Verhältnis heraus und führen Sie es aus. Holen Sie sich das geschätzte Dichteverhältnis zu den Testdaten.
compute_density_ratio = get_method(result, "compute_density_ratio")
w_hat = compute_density_ratio(r_y)
#Dichteverhältnis in Anomalie umwandeln
anom_score = -np.log(w_hat)

Sie können sehen, dass mehr Testdaten den abnormalen Schwellenwert überschreiten.

img003.png

Vergleichen wir die Originaldaten mit dem Grad der Anomalie, zählen die Anzahl der Daten, die den Schwellenwert für den Grad der Anomalie überschreiten, und betrachten die akkumulierten Daten. (Wenn Sie jedoch das auf den Schwellenwert festgelegte Perzentil subtrahieren, wird der Durchschnitt zu normalen Zeiten 0 sein.) Wenn Sie diese Daten als Zeitreihendaten betrachten, ist das erste fast normal und überschreitet den Schwellenwert nicht sehr. Sie können sehen, dass sich die Tendenz der Daten der zweiten Hälfte geändert hat und die Anzahl stark abnormaler Daten plötzlich zunimmt.

img004.png

</ i> Für zweidimensionale Eingabevariablen

Versuchen wir nun den Fall, in dem die Eingabevariable zweidimensional ist. Die Daten stammen von @ oshokawas Folie https://speakerdeck.com/oshokawa/mi-du-bi-tui-ding-niyoruyi-chang-jian-zhi Ich habe es mit den gleichen Spezifikationen wie erstellt.

Python


#Generierung von Trainingsdaten
rd.seed(71)

cov_n = [[ 10,  0],
         [  0, 10]]
M = 2
m_n_1 = (10, 10)
m_n_2 = (20, 20)
n_data = 500

X_norm = np.empty((2*n_data, M))
X_norm[:n_data] = st.multivariate_normal.rvs(mean=m_n_1, cov=cov_n, size=n_data)
X_norm[n_data:] = st.multivariate_normal.rvs(mean=m_n_2, cov=cov_n, size=n_data)
N_norm = 2*n_data

#Generierung von Testdaten
cov_n_1 = [[ 5,  0],
           [  0, 5]]
cov_n_2 = [[ .1,  0],
           [  0,.1]]

m_n_1 = (12.5, 12.5)
m_n_2 = (17.5, 17.5)
m_n_3 = (15, 15)
m_n_4 = (17.5, 12.5)
m_n_5 = (12.5, 17.5)

n_data_a = 50
N_anom = 5*n_data_a
X_anomaly = np.empty((n_data_a*5, M))
for i, m, c in zip(range(5), [m_n_1,m_n_2,m_n_3,m_n_4,m_n_5], [cov_n_1,cov_n_1,cov_n_2,cov_n_1,cov_n_1]):
    X_anomaly[n_data_a*i:n_data_a*(i+1)] = st.multivariate_normal.rvs(mean=m, cov=c, size=n_data_a)

X = np.r_[X_norm, X_anomaly]

Die Daten sehen so aus. Eine Maschine hat einen Rotationsmodus mit niedriger Geschwindigkeit und einen Rotationsmodus mit hoher Geschwindigkeit, und der Zustand ändert sich in der Mitte, aber es wird mit dem Bild von Daten gemacht, dass die Bewegung gegen Ende verdächtig ist.

img005.png

img007.png

Wie bei eindimensionalen Daten wird es in einen R-Datenrahmen konvertiert und dann einem Dens-Verhältnis unterzogen.

Python


#In R-Datenrahmen konvertieren
r_x = pandas2ri.py2ri(pd.DataFrame(X_norm))
r_y = pandas2ri.py2ri(pd.DataFrame(X_anomaly))
r_xall = pandas2ri.py2ri(pd.DataFrame(X))

#Führen Sie eine Schätzung des Dichteverhältnisses durch
result = densratio.densratio(r_x, r_y)

out


//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: ################## Start uLSIF ##################

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: Searching optimal sigma and lambda...

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.001, lambda = 0.001, score = 0.000

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.001, lambda = 3.162, score = -0.000

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.032, lambda = 0.001, score = -0.000

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.100, lambda = 0.001, score = -0.002

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 0.316, lambda = 0.001, score = -0.352

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning:   sigma = 1.000, lambda = 0.010, score = -2.466

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: Found optimal sigma = 1.000, lambda = 0.010.

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: Optimizing alpha...

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: End.

  warnings.warn(x, RRuntimeWarning)
//anaconda/lib/python3.5/site-packages/rpy2/rinterface/__init__.py:185: RRuntimeWarning: ################## Finished uLSIF ###############

  warnings.warn(x, RRuntimeWarning)

Der Ablauf, das Dichteverhältnis herauszunehmen, den Grad der Abnormalität zu berechnen und den Schwellenwert aus den Trainingsdaten festzulegen, ist auch diesmal der gleiche!

Python


#Holen Sie sich das geschätzte Dichteverhältnis
compute_density_ratio = get_method(result, "compute_density_ratio")
w_hat = np.asanyarray(compute_density_ratio(r_xall))

#Abnormalitätsschwelle
anom_percentile = 5 #Untere 5 des geschätzten Dichteverhältnisses%Das Folgende ist abnormal

#Legen Sie einen Schwellenwert für abnormale Trainingsdaten unterhalb des angegebenen Perzentils fest
anom_score = w_hat[:len(X_norm)]
thresh = np.percentile(anom_score, anom_percentile)

Zeichnen wir das Dichteverhältnis. img008.png

Konvertieren Sie das Dichteverhältnis der Trainingsdaten in den Grad der Anomalie und erstellen Sie ein Diagramm. In ähnlicher Weise wird die Abnormalitätsbeurteilung unter Verwendung des angegebenen Perzentils der Trainingsdaten (diesmal 95% Punkt) als Schwellenwert durchgeführt. img009.png

Zählen wir anhand des berechneten Schwellenwerts den Grad der Abnormalität aller Daten und die Anzahl der Daten, die den Schwellenwert überschreiten, und sehen die Tendenz. Immerhin nimmt der Grad der Abnormalität zu, wenn die Situation in der zweiten Hälfte seltsam ist. img010.png

Schauen wir uns die letzte geschätzte Dichte auf der Wärmekarte an. Das linke ist das geschätzte Dichteverhältnis, das mittlere ist die Dichte der Trainingsdaten und das rechte ist die Dichte der Testdaten. img011.png

</ i> Referenz

Vollständiger Quellcode für diese Seite (Die defekten Teile in diesem Artikel, wie z. B. der Zeichnungscode für Grafiken, werden hier ebenfalls beschrieben.) 　https://github.com/matsuken92/Qiita_Contents/blob/master/General/Density_ratio_R.ipynb

rpy2-Paketreferenz 　http://rpy2.readthedocs.io/

Densratio-Paket (@hoxo_m) 　　https://github.com/hoxo-m/densratio

Erkennung von Abnormalitäten durch Schätzung des Dichteverhältnisses (@oshokawa) 　https://speakerdeck.com/oshokawa/mi-du-bi-tui-ding-niyoruyi-chang-jian-zhi

"Erkennung von Anomalien und Erkennung von Veränderungen" von Tsuyoshi Ide und Masaru Sugiyama (Machine Learning Professional Series) 　https://www.amazon.co.jp/dp/4061529080

Singularitätserkennungsverfahren unter Verwendung einer Dichteverhältnisschätzung 　http://2boy.org/~yuta/publications/ibis2007.pdf

Recommended Posts