Utiliser le package d'estimation du ratio de densité R densratio de Python

Cette fois, l'article est une combinaison des deux thèmes, l'un est d'appeler le package R de Python et de l'utiliser, et l'autre est d'essayer de détecter des anomalies en utilisant l'estimation du rapport de densité R package densratio. est. R est très attractif car il y a tellement de packages de méthodes statistiques publiés (plus de 8000), et certains d'entre eux ne sont pas en Python. Cependant, si vous êtes habitué à Python, vous souhaiterez peut-être prétraiter les données et dessiner des graphiques dans le Python familier. (Il s'agit de moi: sweat_smile :) Pour ces personnes, voici comment utiliser le package R de Python.

Le texte complet du code source Python peut être trouvé à [ici] sur GitHub (https://github.com/matsuken92/Qiita_Contents/blob/master/General/Density_ratio_R.ipynb), veuillez donc l'utiliser comme il convient.

</ i> Environnement

L'environnement que j'ai essayé est le suivant. J'essaye sur Mac et Anaconda. Si Anaconda est installé, aucune préparation particulière n'est requise. Le contenu de cet article est essentiellement vérifié avec Jupyter Notebook (iPython Notebook).

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> Appeler le package R depuis Python

Préparation

Tout d'abord, il faut que R soit installé. Ensuite, démarrez R (ou RStudio) une fois et installez le package dens ratio comme indiqué ci-dessous. Il est utilisé pour l'estimation du rapport de densité. Une fois l'installation terminée, vous pouvez quitter R.

r


install.packages("densratio")

Accéder à R depuis Python

Installez une bibliothèque appelée rpy2 pour accéder à R depuis Python.

bash


pip install rpy2

Nous accèderons à R en utilisant ce rpy2. Tout d'abord, importez diverses choses. Vous pouvez désormais utiliser des objets R et effectuer une conversion entre les dataframes Pandas et les dataframes R.

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> Essayez le package d'estimation du ratio de densité dens ratio

Je vais essayer l'estimation du rapport de densité en utilisant le package densratio créé par boss @hoxo_m et la détection d'anomalies en l'utilisant.

Je vais essayer de charger immédiatement le package de rapport R dens avec Python, mais je peux l'importer simplement en écrivant une ligne comme indiqué ci-dessous.

Python


densratio = importr("densratio")

Essayez avec des données unidimensionnelles

L'estimation du rapport de densité est une méthode d'estimation directe du rapport de densité au lieu d'estimer chaque distribution, puis de prendre le rapport lorsque la distribution des variables d'entrée appelées décalage de covariable est différente entre les données d'apprentissage et les données de test. (Une covariable est une variable d'entrée)

Données d'entraînement $ \ mathcal {D} = \ {\ boldsymbol {x} ^ {(1)}, \ cdots, \ boldsymbol {x} ^ {(N)} \} $ sont des données normales uniquement (si anormales) Même s'il y a des données, il s'agit d'un très petit nombre), et testez les données $ \ mathcal {D} '= \ {\ boldsymbol {x} ^ {' (1)}, \ cdots, \ boldsymbol {x} ^ {' (N ')} \} $ contient un certain nombre de données anormales, et ces données anormales sont détectées à partir du rapport de densité.

En supposant que la distribution des données d'entraînement est $ p_ {tr} (\ boldsymbol {x}) $ et que la distribution des données de test est $ p_ {te} (\ boldsymbol {x}) $, quel est le rapport de densité?

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

C'est une valeur exprimée sous la forme. Estimer $ \ hat {w} (\ boldsymbol {x}) $ pour ce $ w (\ boldsymbol {x}) $

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

Et il est représenté par une combinaison linéaire de la fonction de base $ \ {\ varphi_l (\ boldsymbol {x}) \} _ {l = 1} ^ {N} $. Ici, $ \ varphi_l (\ boldsymbol {x}) $ est appelé une fonction de base, et cette fois le noyau RBF est utilisé pour cette fonction de base.

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

Est utilisé. Puisque $ \ boldsymbol {x} ^ {l} $ représente les $ l $ th données d'entraînement, cela ressemble à une distribution normale avec les constantes de normalisation centrées sur les données d'entraînement omises.

Ce qui suit est un exemple, en supposant qu'il existe 5 données d'entraînement «[0,5, 4,5, 5,0, 7,0, 9,0]», à partir desquelles «α = [.1, .25, .35, .2, .1] Les noyaux pondérés sont additionnés, et la partie inférieure montre comment il est représenté par le rapport de densité estimé $ \ hat {w} (\ boldsymbol {x}) $.

img012.png

Donc, cette estimation du rapport de densité consiste à déterminer le rapport de mélange du noyau $ \ boldsymbol {\ alpha} $. $ \ Sigma $, qui indique la largeur du noyau, est un hyper paramètre, qui est déterminé par validation croisée.

Veuillez consulter les livres, pages et articles mentionnés ci-dessus pour savoir pourquoi cela fonctionne et comment décider de $ \ boldsymbol {\ alpha} $.

Génération de données artificielles

Le problème habituel des valeurs aberrantes consiste à calculer le degré d'anomalie individuellement et à le comparer au seuil, en supposant qu'une nouvelle donnée est obtenue, mais cette méthode est caractérisée par le fait que les données de test sont des échantillons multiples. .. Premièrement, ces données sont générées artificiellement à l'aide de nombres aléatoires.

De plus, puisque densratio est un package R, utilisez pandas2ri.py2ri pour transférer des données du monde Pandas vers le monde R.

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)

#Convertissez les dataframes Python en dataframes R!
r_x = pandas2ri.py2ri(pd.DataFrame(x))
r_y = pandas2ri.py2ri(pd.DataFrame(y))

Les données anormales ont une moyenne légèrement plus élevée et une variance plus grande que les données normales.

img001.png

Effectuer une estimation du rapport de densité

Le rapport de densité est calculé en utilisant le rapport de densité pour ces données. Deux types de méthodes d'estimation du rapport de densité peuvent être utilisés dans ce progiciel, KLIEP, qui estime le rapport de densité à l'aide de la quantité d'informations Calback Libra, et uLSIF, qui utilise la méthode d'estimation du rapport de densité carré minimum sans contrainte. Cette fois, nous utiliserons uLSIF, qui est également la valeur par défaut et peut effectuer des calculs à grande vitesse. Insérez simplement le r_x, r_y converti en trame de données R dans cette fonction. C'est pratique: sourire: (Cette fois, $ \ sigma $ a été mis à 0,2 sans validation croisée)

Python


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

Un message d'avertissement est affiché pendant ce processus, mais il s'agit d'un journal du processus de détermination des hyper paramètres par validation croisée au lieu d'avertissement, mais il est affiché comme avertissement dans les spécifications rpy2.

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)

Extrayez et exécutez compute_density_ratio, qui est une fonction pour extraire le rapport de densité du résultat, côté Python. Tout d'abord, obtenez le rapport de densité estimé w_hat à l'échantillon normal. De plus, si le rapport de densité est $ \ hat {w} $, le degré d'anomalie est moins le logarithme.

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

Puisqu'elle est exprimée par, cette conversion est effectuée et calculée.

De plus, je voudrais faire un jugement d'anomalie en utilisant le 99% percentile de la distribution de ces données d'entraînement comme seuil d'anomalie.

Python


#Calculer et extraire le rapport de densité
compute_density_ratio = get_method(result, "compute_density_ratio")
w_hat = compute_density_ratio(r_x)

#Convertir le rapport de densité en anomalie
anom_score = -np.log(w_hat)
thresh = np.percentile(anom_score, 99)    #99 du rapport de densité des échantillons normaux%Définissez le percentile comme seuil pour juger les anomalies

Voici un graphique de ces «w_hat» et «anom_score».

img002.png

Ensuite, calculez et dessinez le rapport de densité des données de test et l'histogramme du degré d'anomalie.

Python


#Calculer à partir du résultat_density_Retirez le ratio et exécutez-le. Obtenez le rapport de densité estimé aux données de test.
compute_density_ratio = get_method(result, "compute_density_ratio")
w_hat = compute_density_ratio(r_y)
#Convertir le rapport de densité en anomalie
anom_score = -np.log(w_hat)

Vous pouvez voir qu'il y a plus de données de test qui dépassent le seuil anormal.

img003.png

Comparons les données d'origine avec le degré d'anomalie, comptons le nombre de données qui dépasse le degré de seuil d'anomalie et examinons les données accumulées. (Cependant, en soustrayant le percentile défini au seuil, la moyenne sera de 0 en temps normal.) Si vous considérez ces données comme des données de série chronologique, la première est presque normale et ne dépasse pas beaucoup le seuil. , Vous pouvez voir que la tendance de la dernière moitié des données a changé et que le nombre de données hautement anormales augmente soudainement.

img004.png

</ i> Pour les variables d'entrée bidimensionnelles

Essayons maintenant le cas où la variable d'entrée est bidimensionnelle. Les données proviennent de la diapositive de @ oshokawa https://speakerdeck.com/oshokawa/mi-du-bi-tui-ding-niyoruyi-chang-jian-zhi Je l'ai créé avec les mêmes spécifications que.

Python


#Génération de données d'entraînement
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

#Génération de données de test
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]

Les données ressemblent à ceci. Une machine a un mode de rotation à basse vitesse et un mode de rotation à grande vitesse, et l'état change au milieu, mais il est fait avec l'image des données que le mouvement est suspect vers la fin.

img005.png

img007.png

Comme pour les données unidimensionnelles, elles sont converties en une trame de données R puis soumises au rapport de densité.

Python


#Convertir en trame de données R
r_x = pandas2ri.py2ri(pd.DataFrame(X_norm))
r_y = pandas2ri.py2ri(pd.DataFrame(X_anomaly))
r_xall = pandas2ri.py2ri(pd.DataFrame(X))

#Effectuer une estimation du rapport de densité
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)

Le déroulement de la prise du rapport de densité, du calcul du degré d'anomalie et du réglage du seuil à partir des données d'entraînement est également le même cette fois!

Python


#Obtenez le rapport de densité estimé
compute_density_ratio = get_method(result, "compute_density_ratio")
w_hat = np.asanyarray(compute_density_ratio(r_xall))

#Seuil d'anomalie
anom_percentile = 5 #Côté inférieur du rapport de densité estimé 5%Ce qui suit est anormal

#Définir un seuil pour les données d'entraînement anormales en dessous du centile spécifié
anom_score = w_hat[:len(X_norm)]
thresh = np.percentile(anom_score, anom_percentile)

Tracons le rapport de densité. img008.png

Convertissez le rapport de densité des données d'entraînement en degré d'anomalie et créez un graphique. De même, le jugement d'anomalie est effectué en utilisant le percentile spécifié des données d'apprentissage (95% cette fois) comme seuil. img009.png

Sur la base du seuil calculé, comptons le degré d'anomalie de toutes les données et le nombre de données qui dépasse le seuil et voyons la tendance. Après tout, le degré d'anomalie augmente là où la situation dans la seconde moitié est étrange. img010.png

Regardons la dernière densité estimée sur la carte thermique. La gauche est le rapport de densité estimé, le milieu la densité des données d'entraînement et la droite la densité des données de test. img011.png

</ i> Référence

Code source complet pour cette page (Les parties cassées de cet article, telles que le code de dessin graphique, sont également décrites ici.) 　https://github.com/matsuken92/Qiita_Contents/blob/master/General/Density_ratio_R.ipynb

référence du package rpy2 　http://rpy2.readthedocs.io/

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

Détection d'anomalies par estimation du rapport de densité (@oshokawa) 　https://speakerdeck.com/oshokawa/mi-du-bi-tui-ding-niyoruyi-chang-jian-zhi

"Détection d'anomalies et détection de changement" par Tsuyoshi Ide et Masaru Sugiyama (série professionnelle d'apprentissage automatique) 　https://www.amazon.co.jp/dp/4061529080

Méthode de détection de singularité utilisant l'estimation du rapport de densité 　http://2boy.org/~yuta/publications/ibis2007.pdf

Recommended Posts