[PYTHON] Introduction à la détection des anomalies et résumé des méthodes

――Un comportement anormal dépend de problèmes individuels, tels qu'une augmentation ou une diminution drastique des valeurs observées dans la même section ou un changement significativement plus petit. --La principale différence par rapport aux anomalies statiques est que les anomalies dynamiques n'ont pas toujours des valeurs anormales. Par exemple, la partie centrale du graphique ci-dessous semble anormale, mais la valeur elle-même est une valeur possible. Dans ce cas, vous trouverez des anomalies dans le comportement des «changements de fréquence».

import numpy as np
import matplotlib.pyplot as plt

x1 = np.concatenate([np.linspace(-10, -3, 50), np.zeros(50), np.linspace(3, 9, 50)])
x2 = np.concatenate([np.zeros(50), np.linspace(-3, 3, 50), np.zeros(50)])
y = np.sin(x1) + np.sin(5*x2)
plt.plot(np.linspace(-10, 9, 150), y)

Exprimé sous forme de formule, chaque point $ x_j $ est $ p (x) = \ int p (x, t) {\ rm d} t $ dans les données de la série $ \ {x_1, ..., x_t \} $ Même si elle est normale au sens de, la probabilité conditionnelle du temps $ p (x_j \ mid t) = \ int \ cdots \ int p (x_1, ... x_t \ mid t) {\ rm d} t_1 \ cdots { \ rm d} t_ {j-1} {\ rm d} t_ {j + 1} \ cdots {\ rm d} t_t $ peut être anormal.

Cependant, même dans ce cas, la partie centrale peut ne pas être anormale. Par exemple, il est possible que le graphique supérieur soit le signal de fonctionnement de la machine de l'usine et que le mode de fonctionnement soit commuté uniquement dans la partie centrale, ce n'est donc pas anormal.

Comme mentionné ci-dessus, il existe différents aspects des anomalies, mais la référence [1] résume les types d'anomalies de cette manière:

出典:[1]

L'anomalie de point dans la figure en haut à gauche ressemble à des données de bruit d'un seul coup (valeurs aberrantes) dues à des erreurs d'observation, etc., mais l'anomalie de groupe ressemble plus à des groupes de données générés à partir d'une distribution différente qu'à des valeurs aberrantes. On suppose que les deux nécessitent des approches de détection d'anomalies différentes. Lorsque vous atteignez Anomalie sémantique en bas à droite, il est difficile de la gérer avec la méthode classique de détection des anomalies, et une approche d'apprentissage en profondeur est nécessaire.

Comme vous pouvez le voir, il n'y a pas de solution unifiée pour «anomalie», et «anomalie» doit être définie au cas par cas.

L'algorithme de détection des anomalies est essentiellement un apprentissage non supervisé

L'algorithme de détection des anomalies utilise essentiellement un apprentissage non supervisé plutôt qu'un apprentissage supervisé. Il y a deux raisons:

1. Principalement parce qu'il y a peu de données anormales
2. Les données anormales sont diverses et difficiles à modéliser de manière exhaustive

est. Des centaines à des milliers de données d'anomalies ne doivent pas être collectées, en particulier dans les sites où les anomalies sont critiques (médicales ou en usine). De plus, la raison 2 est celle expliquée au début. Il est plus facile et plus simple de deviner «normal / pas normal» que de deviner «normal / anormal». Ce dernier ne mentionne pas de données anormales et ne peut être obtenu qu'en comprenant la structure des données normales. L'apprentissage non supervisé est un cadre permettant de comprendre la structure des données en estimant $ p (x) $ à partir de l'échantillon $ x_i (i = 1, ..., n) $, donc la détection des anomalies est fondamentalement non supervisée. Vous utiliserez l'apprentissage. Cependant, il doit être possible de supposer que l'échantillon pour l'apprentissage non supervisé ne contient pas de données anormales ou est extrêmement petit.

En passant, la référence [2] indique ce qui suit et est ma préférée:

Un passage du roman de Tolstoï "Anna Karénine" qui dit: "Les familles heureuses se ressemblent toutes; chaque famille malheureuse est malheureuse à sa manière" Est suggestif.

Les algorithmes de détection d'anomalies sont à peu près divisés en quatre types

Ensuite, je présenterai concrètement l'algorithme de détection des anomalies.

Les algorithmes peuvent être globalement classés en quatre approches. Approche de classification, approche probabiliste, approche de reconstruction, approche à distance, respectivement. En plus de ces quatre classifications, voici un tableau qui classe s'il faut ou non utiliser le deep learning:

出典：[1]

Chaque mot tel que PCA ou OC-SVM est le nom de l'algorithme de détection d'anomalies qui a été proposé dans le passé. Les détails de chaque algorithme seront présentés plus tard, nous ne donnerons donc ici qu'un aperçu de chacune des quatre approches. (Je ne comprends pas entièrement le tableau ci-dessous, je vous serais donc reconnaissant si vous pouviez me donner votre avis.)

Le tableau ci-dessus a été créé en référence aux matériaux suivants.

1. Hido: Introduction à Jubatus Casual Talks # 2 Détection d'anomalies
2. PANG, et al. Deep Learning for Anomaly Detection: A Review

** §2 Introduction et expérimentation de l'algorithme **

Résumé des méthodes à introduire

• Classification Approach
  • One-Class SVM
• Probability Approach --Méthode T2 de l'hôtellerie
  • GMM
  • KDE
• Reconstruction Approach
  • AE
• Distance Approach
  • k-NN
  • LOF

Classification Approach One-Class SVM Supposons que vous receviez les données $ D = \ {x_1, ..., x_N \} $. Cependant, on suppose que $ D $ ne contient pas de données anormales, ou si c'est le cas, il est extrêmement petit.

Dans le chapitre précédent, nous avons introduit que l'approche de classification estime la frontière qui sépare normal et anormal. ** SVM à une classe est un algorithme qui calcule et borde une sphère qui entoure complètement les données normales. ** (Une sphère n'est pas nécessairement tridimensionnelle.)

Si vous rendez le rayon très grand, vous pouvez créer une sphère qui contient complètement $ D $, mais si la sphère est trop grande, False Negative sera 1 et je ne suis pas content. Donc,

-Rendre le volume de la sphère aussi petit que possible -Autoriser certaines données en $ D $ à s'étendre au-delà de la sphère

Calculez une sphère à usage général en autorisant deux points. Plus précisément, le rayon optimal de la sphère $ R $ et le centre $ b $ peuvent être obtenus en résolvant le problème d'optimisation suivant:

python

\min_{R,b,u} \left\{R^2 + C\sum_{n=1}^{N}u_n \right\} \quad {\rm s.t.} \quad \|x_n - b\|^2 \le R^2 + u_n

Cependant, $ u $ est un terme de pénalité qui s'étend au-delà de la sphère, et $ C $ est un paramètre de l'importance de cette pénalité. La conversion non linéaire des données par les fonctions du noyau est possible avec SVM, mais bien sûr, elle est également disponible avec SVM à une classe. Une bonne utilisation des fonctions du noyau vous donnera la bonne sphère.

Faisons une expérience numérique. Les données ont été préparées comme suit.

python

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

sns.set()
np.random.seed(123)

#Moyenne et variance
mean1 = np.array([2, 2])
mean2 = np.array([-2,-2])
cov = np.array([[1, 0], [0, 1]])

#Génération de données normales(Généré à partir de deux distributions normales)
norm1 = np.random.multivariate_normal(mean1, cov, size=50)
norm2 = np.random.multivariate_normal(mean2, cov, size=50)
norm = np.vstack([norm1, norm2])

#Génération de données anormales(Généré à partir d'une distribution uniforme)
lower, upper = -10, 10
anom = (upper - lower)*np.random.rand(20,20) + lower

s=8
plt.scatter(norm[:,0], norm[:,1], s=s, color='green', label='normal')
plt.scatter(anom[:,0], anom[:,1], s=s, color='red', label='anomaly')
plt.legend()

python

from sklearn import svm
xx, yy = np.meshgrid(np.linspace(-10, 10, 500), np.linspace(-10, 10, 500))

# fit the model
clf = svm.OneClassSVM(nu=0.1, kernel="rbf", gamma=0.1)
clf.fit(norm)

def draw_boundary(norm, anom, clf):
    # plot the line, the points, and the nearest vectors to the plane
    Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)

    plt.title("One-Class SVM")
    plt.contourf(xx, yy, Z, levels=np.linspace(Z.min(), Z.max(), 7), cmap=plt.cm.PuBu)
    bd = plt.contour(xx, yy, Z, levels=[0], linewidths=2, colors='darkred')
    plt.contourf(xx, yy, Z, levels=[0, Z.max()], colors='palevioletred')

    nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
    am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

    plt.axis('tight')
    plt.xlim((-10, 10))
    plt.ylim((-10, 10))
    plt.legend([bd.collections[0], nm, am],
               ["learned frontier", "normal", "anomaly"],
               loc="upper left",
               prop=matplotlib.font_manager.FontProperties(size=11))
    plt.show()

draw_boundary(norm, anom, clf)

Au fait, si vous réduisez le paramètre $ C $, c'est-à-dire que la pénalité qui dépasse de la sphère est plus claire, la limite changera comme ceci. (Dans sklearn, il est représenté par le paramètre $ \ nu $, mais il semble que ce soit dans le dénominateur, il semble donc que $ \ nu $ est rendu plus petit et $ = C $ est rendu plus grand.)

python

# fit the model
clf = svm.OneClassSVM(nu=0.5, kernel="rbf", gamma=0.1)
clf.fit(norm)
draw_boundary(norm, anom, clf)

D'autres propositions incluent One-Class NN. Dans One-Class SVM, je cherchais une transformation non linéaire appropriée en ajustant la fonction du noyau, mais l'idée est de la remplacer et de l'automatiser avec NN. L'implémentation semble ennuyeuse, donc je ne l'ai pas fait.

Probability Approach

Méthode hôtelière T2

La méthode T2 d'Hoteling est une méthode classique et a été théoriquement analysée. En conclusion, ** En supposant une distribution normale, cela devient un algorithme ** qui utilise la propriété que le score d'anomalie suit la distribution $ \ chi ^ 2 $.

Par souci de simplicité, nous utiliserons ici des données unidimensionnelles. Comme d'habitude, étant donné les données $ D $. Ici, chaque donnée $ x $ est normalement distribuée

python

p(x\mid \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{1}{2\sigma^2}(x-\mu)^2\right\}

On suppose qu'il suit **. En ajustant cette distribution normale aux données par la méthode d'estimation la plus probable

python

\begin{align*}
    \hat{\mu} &= \frac{1}{N}\sum_{n=1}^{N}x_n \\
    \hat{\sigma}^2 &= \frac{1}{N}\sum_{n=1}^{N}(x_n - \hat{\mu})^2
\end{align*}

Et chaque paramètre est estimé. En d'autres termes, nous avons estimé que les données $ D $ suivent la distribution normale $ p (x \ mid \ hat {\ mu}, \ hat {\ sigma} ^ 2) $.

Maintenant, introduisons le score d'anomalie. L'anomalie est une statistique qui indique littéralement le degré d'anomalie, et la probabilité logarithmique négative est souvent utilisée.

python

prob = np.linspace(0.01, 1, 100)
neg_log_likh = [-np.log(p) for p in prob]

plt.plot(prob, neg_log_likh)
plt.title('Negative Log Likelihood')
plt.ylabel('Anomaly Score')
plt.xlabel('Prob.')

Comme le montre le graphique, la vraisemblance logarithmique négative a une forme qui prend une grande valeur là où la probabilité est faible, et elle semble fonctionner comme une anomalie. Si vous calculez l'anomalie de la distribution normale que vous venez d'estimer

python

-\log p(x\mid \hat{\mu}, \hat{\sigma}^2) = \frac{1}{2\hat{\sigma}^2}(x-\hat{\mu})^2 + \frac{1}{2}\log(2\pi\hat{\sigma}^2)

Ce sera. Puisque l'anomalie est une quantité définie pour les données $ x $, l'anomalie $ a (x) $ l'est, à l'exclusion des termes non liés à $ x $.

python

a(x) = \left( \frac{x - \hat{\mu}}{\hat{\sigma}} \right)^2

Est calculé. La molécule d'anomalie signifie que «plus les données sont éloignées de la moyenne, plus l'anomalie est élevée». Le dénominateur exprime "combien de données sont dispersées?", Et le sentiment que les données dispersées sont tolérées même si elles sont éloignées de la moyenne, et inversement, les données denses ne permettent aucun écart. Il y a.

En fait, on sait que cette anomalie $ a (x) $ suit la distribution du chi carré $ \ chi ^ 2_1 $ avec un degré de liberté de $ 1 $ lorsque la taille de l'échantillon $ N $ est suffisamment grande.

python

from scipy.stats import chi2
x = np.linspace(0, 8, 1000)
for df in [1,2,3,4,5]:
    plt.plot(x, chi2.pdf(x ,df = df), label='deg. of freedom:'+str(df))
plt.legend()
plt.ylim(0, 0.5)
plt.title('Chi-Square Dist.')

Après cela, si vous décidez du seuil $ a_ {th} $ au niveau de signification de $ 5 $%, vous pouvez détecter des anomalies. Sommaire,

1. Supposons que les données suivent une distribution normale.
2. Estimer les paramètres de distribution normale $ \ mu $, $ \ sigma $ par la méthode d'estimation la plus probable
3. À partir de la distribution du chi carré, déterminez le point $ a_ {th} $ qui sépare normal et anormal.
4. Si les nouvelles données $ x_ {new} $ sont $ a (x_ {new})> a_ {th} $, elles sont considérées comme anormales.

Est la méthode T2 de l'hôtellerie. Cette fois, je l'ai expliqué dans une dimension, mais le degré d'anomalie suit également la distribution du chi carré en plusieurs dimensions. Cependant, le degré de liberté = le nombre de dimensions.

Expérimentons. Les données ont été créées comme suit.

python

norm = np.random.normal(0, 1, 20) #Données normales
anom = np.array([-3, 3]) #Données anormales
plt.scatter(norm, [0]*len(norm), color='green', label='normal')
plt.scatter(anom, [0]*len(anom), color='red', label='anomaly')

python

#Estimation la plus probable
N = len(norm)
mu_hat = sum(norm)/N
s_hat = sum([(x-mu_hat)**2 for x in norm])/N

#Utilisez le point de division de la distribution du chi carré comme seuil
a_th_5 = chi2.ppf(q=0.95, df=1)
a_th_30 = chi2.ppf(q=0.7, df=1)

#Calculez le degré d'anomalie de chaque donnée
data = np.hstack([norm, anom])
anom_score = []
for d in data:
    score = ((d - mu_hat)/s_hat)**2
    anom_score.append(score)

#dessin
norm_pred_5 = [d for d,a in zip(data, anom_score) if a<=a_th_5]
anom_pred_5 = [d for d,a in zip(data, anom_score) if a>a_th_5]
norm_pred_30 = [d for d,a in zip(data, anom_score) if a<=a_th_30]
anom_pred_30 = [d for d,a in zip(data, anom_score) if a>a_th_30]

fig, axes = plt.subplots(1,3, figsize=(15,5))

axes[0].scatter(norm, [0]*len(norm), color='green', label='normal')
axes[0].scatter(anom, [0]*len(anom), color='red', label='anomaly')
axes[0].set_title('Truth')
axes[1].scatter(norm_pred_5, [0]*len(norm_pred_5), color='green', label='normal pred')
axes[1].scatter(anom_pred_5, [0]*len(anom_pred_5), color='red', label='anomaly pred')
axes[1].set_title('Threshold:5%')
axes[2].scatter(norm_pred_30, [0]*len(norm_pred_30), color='green', label='normal pred')
axes[2].scatter(anom_pred_30, [0]*len(anom_pred_30), color='red', label='anomaly pred')
axes[2].set_title('Threshold:30%')
plt.legend()

On constate que si le seuil d'anomalie est réduit, la zone considérée comme anormale s'agrandira. Le seuil doit être déterminé en fonction du problème. Par exemple, le seuil doit être défini petit lorsque la négligence des anomalies est critique, comme dans le domaine médical, et peut être fixé un peu plus haut dans le cas du nettoyage des données.

GMM(Gaussian Mixture Model) La méthode T2 d'Hoteling a supposé une distribution normale dans les données, mais dans certains cas, cela peut être inapproprié. Par exemple, si les données sont concentrées à plusieurs endroits (plusieurs pics au lieu de pics uniques) comme indiqué ci-dessous, ajustement avec une distribution normale

python

lower, upper = -3, -1
data1 = (upper - lower)*np.random.rand(20) + lower
lower, upper = 1, 3
data2 = (upper - lower)*np.random.rand(30) + lower
data = np.hstack([data1, data2])

#Estimation la plus probable
N = len(data)
mu_hat = sum(data)/N
s_hat = sum([(x-mu_hat)**2 for x in data])/N

#terrain
from scipy.stats import norm
X = np.linspace(-7, 7, 1000)
Y = norm.pdf(X, mu_hat, s_hat)

plt.scatter(data, [0]*len(data))
plt.plot(X, Y)

On estime que le pic de la distribution normale arrive à l'endroit où aucune donnée n'est observée. Dans de tels cas, il est préférable de définir une distribution légèrement plus compliquée au lieu de la distribution normale. Cette fois, nous allons introduire la détection d'anomalies ** par GMM (Gaussian Mixture Model) à titre d'exemple.

GMM est une distribution qui exprime un type multi-pic en ajoutant plusieurs distributions normales. Le nombre à ajouter est souvent exprimé en $ K $ et est appelé le nombre de composants. Les formules et les visuels GMM sont les suivants (une dimension).

python

p(x\mid \pi, \mu, \sigma) = \sum_{k=1}^{K}\pi_k N(x\mid \mu_k, \sigma_k)

cependant,

python

0 \le\pi_k\le 1,\quad \sum_{k=1}^{K}\pi_k=1,\quad N(x\mid\mu,\sigma)={\distribution normale du texte}

est.

python

#terrain
from scipy.stats import norm
X = np.linspace(-7, 7, 1000)

Y1 = norm.pdf(X, 1, 1)
Y2 = norm.pdf(X, -3, 2)
Y = 0.5*Y1 + 0.5*Y2 

plt.plot(X, Y1, label='Component1')
plt.plot(X, Y2, label='Component2')
plt.plot(X, Y,  label='GMM')
plt.legend()

La moyenne ($ \ mu ) / variance ( \ sigma $) et le rapport de mélange $ \ pi $ de chaque distribution normale sont des paramètres appris à partir des données. Cependant, le nombre de composants $ K $ doit être fixé. Les paramètres que l'analyste décide de cette manière sont appelés hyper paramètres.

Expérimentons la détection d'anomalies par GMM. Observons également comment la forme de la distribution estimée change en fonction du nombre de composants $ K $. Les données utilisent le même que le SVM à une classe:

python

def draw_GMM(norm, anom, clf, K):
    # plot the line, the points, and the nearest vectors to the plane
    xx, yy = np.meshgrid(np.linspace(-10, 10, 1000), np.linspace(-10, 10, 1000))
    Z = np.exp(clf.score_samples(np.c_[xx.ravel(), yy.ravel()]))
    Z = Z.reshape(xx.shape)

    plt.title("GMM:K="+str(K))
    plt.contourf(xx, yy, Z, levels=np.linspace(Z.min(), Z.max(), 7), cmap=plt.cm.PuBu)
    
    nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
    am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

    plt.axis('tight')
    plt.xlim((-10, 10))
    plt.ylim((-10, 10))
    plt.legend([nm, am],
               ["normal", "anomaly"],
               loc="upper left",
               prop=matplotlib.font_manager.FontProperties(size=11))
    plt.show()

python

#Montage par GMM
from sklearn.mixture import GaussianMixture

for K in [1,2,3,4]:
    clf = GaussianMixture(n_components=K, covariance_type='full')
    clf.fit(norm)
    draw_GMM(norm, anom, clf, K)

En augmentant le nombre de composants de cette manière, vous pouvez voir que même les tendances dans les détails des données sont capturées. Cependant, il faut être prudent car il ne capture pas les détails = les performances de généralisation sont élevées. Le nombre de composants peut être déterminé à l'aide de critères de quantité d'informations tels que WAIC (j'ai essayé d'utiliser WAIC pour le moment: GitHub)

KDE(Kernel Density Estimation) Jusqu'à présent, nous avons supposé une distribution normale et une distribution normale mixte, et estimé les paramètres internes par ajustement aux données. Une telle approche est appelée une approche paramétrique. Mais que faire si la distribution paramétrique ne fonctionne pas?

Une approche non paramétrique est efficace dans de tels cas, et ** KDE est une méthode d'estimation non paramétrique typique **. KDE est une méthode d'estimation de la distribution de probabilité en superposant un noyau (généralement gaussien) sur chaque point de données **. La visualisation est la suivante.

python

from scipy.stats import norm
np.random.seed(123)
lower, upper = -1, 5
data = (upper - lower)*np.random.rand(3) + lower
for i, d in enumerate(data):
    X = np.linspace(lower,upper,500)
    Y = norm.pdf(X, loc=d, scale=1)
    knl = plt.plot(X,Y,color='green', label='kernel'+str(i+1))

# KDE
n = 3
h = 0.24
Y_pdf = []
for x in X:
    Y = 0
    for d in data:
        Y += norm.pdf((x-d)/h, loc=0, scale=1)
    Y /= n*h
    Y_pdf.append(Y)
kde = plt.plot(X,Y_pdf, label='KDE')
plt.scatter(data, [0]*len(data), marker='x')
plt.ylim([-0.05, 1])
plt.legend()

De cette manière, la distribution de probabilité est estimée en superposant les noyaux. Lors de l'utilisation d'une distribution normale pour le noyau, il est nécessaire de définir la distribution (appelée bande passante). Dans les expériences suivantes, observons également le changement de la distribution estimée lorsque cette bande passante est modifiée.

Les détails de la théorie de KDE et de ses forces et faiblesses sont [Suzuki. Data Analysis 10e "Méthode d'estimation de densité non paramétrique"](http://ibis.tu-tokyo.ac.jp/suzuki/lecture/2015/dataanalysis/L9. pdf) était facile à comprendre:

** Distribution estimée **

python

p(x) = \frac{1}{Nh}\sum_{n=1}^{N}K\left(\frac{x-X_n}{h}\right)

Cependant, $ h> 0 $ est la bande passante, $ K () $ est la fonction du noyau et $ X_n $ est chaque point de données. ** Avantages désavantages ** Source: Suzuki. Analyse des données 10e "Méthode d'estimation non paramétrique de la densité"

Expérimentons la détection d'anomalies par KDE. Les données sont comme d'habitude Utilisez le. La distribution estimée lorsque la bande passante est modifiée est la suivante.

python

def draw_KDE(norm, anom, clf, K):
    # plot the line, the points, and the nearest vectors to the plane
    xx, yy = np.meshgrid(np.linspace(-10, 10, 1000), np.linspace(-10, 10, 1000))
    Z = np.exp(clf.score_samples(np.c_[xx.ravel(), yy.ravel()]))
    Z = Z.reshape(xx.shape)

    plt.title("KDE:BandWidth="+str(bw))
    plt.contourf(xx, yy, Z, levels=np.linspace(Z.min(), Z.max(), 7), cmap=plt.cm.PuBu)
    
    nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
    am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

    plt.axis('tight')
    plt.xlim((-10, 10))
    plt.ylim((-10, 10))
    plt.legend([nm, am],
               ["normal", "anomaly"],
               loc="upper left",
               prop=matplotlib.font_manager.FontProperties(size=11))
    plt.show()

python

from sklearn.neighbors import KernelDensity
for bw in [0.2,0.5,1.0,2.0]:
    clf = KernelDensity(bandwidth=bw, kernel='gaussian')
    clf.fit(norm)
    draw_KDE(norm, anom, clf, bw)

Comme vous pouvez le voir, la distribution estimée devient plus fluide à mesure que la bande passante augmente. De plus, étant donné qu'il existe d'autres fonctions du noyau en plus de la distribution normale, il est nécessaire de définir à la fois la "bande passante et la fonction du noyau" de manière appropriée.

Reconstruction Approach Étant donné que l'approche de reconstruction est souvent utilisée pour les données d'image, nous utiliserons également les données d'image ici. Soit les données normales MNIST et les données anormales Fashion-MNIST.

python

#Charge MNIST
import torch
from torch.utils.data import DataLoader
from torchvision import datasets, transforms

#Lecture des données d'entraînement
train_mnist = datasets.MNIST(
    root="./data", 
    train=True, 
    transform=transforms.Compose([
               transforms.ToTensor(),
               transforms.Normalize((0.5,), (0.5,)) #section[-1,1]À l'intérieur, écart type= 0.5, 0.Normalisé à 5
               ]), 
    download=True
)

train_mnist_loader = torch.utils.data.DataLoader(
    train_mnist, batch_size=64, shuffle=True
)

#Lire les données de test
test_mnist = datasets.MNIST(
    root="./data", 
    train=False, 
    transform=transforms.Compose([
               transforms.ToTensor(),
               transforms.Normalize((0.5,), (0.5,)) #section[-1,1]À l'intérieur, écart type= 0.5, 0.Normalisé à 5
               ]), 
    download=True
)

test_mnist_loader = torch.utils.data.DataLoader(
    test_mnist, batch_size=64, shuffle=True
)

python

# Fashion-Charge MNIST
#Lecture des données d'entraînement
train_fashion = datasets.FashionMNIST(
    root="./data", 
    train=True, 
    transform=transforms.Compose([
               transforms.ToTensor(),
               transforms.Normalize((0.5,), (0.5,)) #section[-1,1]À l'intérieur, écart type= 0.5, 0.Normalisé à 5
               ]), 
    download=True
)

train_fashion_loader = torch.utils.data.DataLoader(
    train_fashion, batch_size=64, shuffle=True
)

#Lire les données de test
test_fashion = datasets.FashionMNIST(
    root="./data", 
    train=False, 
    transform=transforms.Compose([
               transforms.ToTensor(),
               transforms.Normalize((0.5,), (0.5,)) #section[-1,1]À l'intérieur, écart type= 0.5, 0.Normalisé à 5
               ]), 
    download=True
)

test_fashion_loader = torch.utils.data.DataLoader(
    test_fashion, batch_size=64, shuffle=True
)

python

import matplotlib.pyplot as plt

#Une fonction qui visualise les images côte à côte à la fois
def imshow(image_set, nrows=2, ncols=10):
    plot_num = nrows * ncols
    fig, ax = plt.subplots(ncols=ncols, nrows=nrows, figsize=(15, 15*nrows/ncols))
    ax = ax.ravel()
    for i in range(plot_num):
        ax[i].imshow(image_set[i].reshape(28, 28), cmap='gray')
        ax[i].set_xticks([])
        ax[i].set_yticks([])

python

#Visualisation d'échantillons MNIST
imshow(train_mnist.data)

python

# Fashion-Visualisation d'échantillons MNIST
imshow(train_fashion.data)

AE(Auto Encoder) AE (Auto Encoder) est un type de NN (Neural Network), et est un NN qui apprend de sorte que l'entrée et la sortie soient égales. Puisqu'il est ennuyeux avec juste une fonction constante, il prend généralement en sandwich une couche intermédiaire avec une dimension plus petite que la dimension d'entrée.

Ce que cela fait, c'est "résumer" les données d'entrée dans la couche intermédiaire. Mathématiquement, cela signifie que les données d'entrée sont mappées sur un espace de faible dimension tout en conservant suffisamment d'informations pour les reconstruire.

** La détection d'anomalies à l'aide d'AE est basée sur l'idée que "les AE formés pour reconstruire des données normales peuvent ne pas être capables de reconstruire correctement les données anormales." ** **

Maintenant, expérimentons la détection d'anomalies par AE. Après la formation AE avec MNIST, nous tenterons de détecter la précision de reconstruction de MNIST et Fashion-MNIST comme le degré d'anomalie.

python

import torch.nn as nn
import torch.nn.functional as F

#Définition de l'EI
class AE(nn.Module):
    def __init__(self):
        super(AE, self).__init__()
        self.fc1 = nn.Linear(28*28, 256)
        self.fc2 = nn.Linear(256, 28*28)
        
    def forward(self, x):
        x = torch.tanh(self.fc1(x))
        x = torch.tanh(self.fc2(x))
        return x
ae = AE()

#Définition de la fonction de perte et méthode d'optimisation
import torch.optim as optim

criterion = nn.MSELoss()
optimizer = optim.Adam(ae.parameters(), lr=0.0001)

python

#Boucle d'apprentissage
N_train = train_mnist.data.shape[0]
N_test = test_mnist.data.shape[0]

train_loss = []
test_loss = []
for epoch in range(10):  
    # Train step
    running_train_loss = 0.0
    for data in train_mnist_loader:
        inputs, _ = data
        inputs = inputs.view(-1, 28*28)

        # zero the parameter gradients
        optimizer.zero_grad()

        # forward + backward + optimize
        outputs = ae(inputs)
        loss = criterion(outputs, inputs)
        loss.backward()
        optimizer.step()

        running_train_loss += loss.item()
    running_train_loss /= N_train
    train_loss.append(running_train_loss)
    
    # Test step
    with torch.no_grad():
        running_test_loss = 0
        for data in test_mnist_loader:
            inputs, _ = data
            inputs = inputs.view(-1, 28*28)
            outputs = ae(inputs)
            running_test_loss += criterion(outputs, inputs).item()
        running_test_loss /= N_test
        test_loss.append(running_test_loss)
    if (epoch+1)%1==0:    # print every 1 epoch
        print('epoch:{}, train_loss:{}, test_loss:{}'.format(epoch + 1, running_train_loss, running_test_loss))
        
plt.plot(train_loss, label='Train loss')
plt.plot(test_loss, label='Test loss')
plt.ylabel('Loss')
plt.xlabel('epoch')
plt.legend()

python

#Visualisons l'état de la reconstruction
with torch.no_grad():
    for data in test_mnist_loader:
        inputs, _ = data
        inputs = inputs.view(-1, 28*28)
        outputs = ae(inputs)

n = 5 # num to viz
fig, axes = plt.subplots(ncols=n, nrows=2, figsize=(15, 15*2/n))
for i in range(n):
    axes[0, i].imshow(inputs[i].reshape(28, 28), cmap='gray')
    axes[1, i].imshow(outputs[i].reshape(28, 28), cmap='gray')
    axes[0, i].set_title('Input')
    axes[1, i].set_title('Reconstruct')
    axes[0, i].set_xticks([])
    axes[0, i].set_yticks([])
    axes[1, i].set_xticks([])
    axes[1, i].set_yticks([])

Il semble qu'il ait été parfaitement reconstruit.

Maintenant, dessinons la distribution de la précision de reconstruction (erreur carrée cette fois) de MNIST et Fashion-MNIST. Si la détection d'anomalie par AE réussit, la précision de Fashion-MNIST doit être faible.

python

#Calcul de la précision de reconstruction
mnist_loss = []
with torch.no_grad():
    for data in test_mnist_loader:
        inputs, _ = data
        for data_i in inputs:
            data_i = data_i.view(-1, 28*28)
            outputs = ae(data_i)
            loss = criterion(outputs, data_i)
            mnist_loss.append(loss.item())

fashion_loss = []
with torch.no_grad():
    for data in test_fashion_loader:
        inputs, _ = data
        for data_i in inputs:
            data_i = data_i.view(-1, 28*28)
            outputs = ae(data_i)
            loss = criterion(outputs, data_i)
            fashion_loss.append(loss.item())

python

plt.hist([mnist_loss, fashion_loss], bins=25, label=['MNIST', 'Fashion-MNIST'])
plt.xlabel('Loss')
plt.ylabel('data count')
plt.legend()
plt.show()

Certes, la reconstruction de Fashion-MNIST ne semble pas bien fonctionner, il semble donc que la détection des anomalies par la précision de reconstruction sera établie.

Enfin, je vais ramasser et dessiner ceux avec une faible perte (= ceux qui sont difficiles à juger comme anormaux) et ceux avec une perte élevée (= ceux qui peuvent être facilement jugés comme anormaux) parmi Fashion-MNIST.

python

lower_loss, upper_loss = 0.1, 0.8
data_low = [] #Stocke les données avec une faible perte
data_up = []  #Stocke les données avec une perte élevée

with torch.no_grad():
    for data in test_fashion_loader:
        inputs, _ = data
        for data_i in inputs:
            data_i = data_i.view(-1, 28*28)
            outputs = ae(data_i)
            loss = criterion(outputs, data_i).item()
            if loss<lower_loss:
                data_low.append([data_i, outputs, loss])
            if loss>upper_loss:
                data_up.append([data_i, outputs, loss])

python

#Index pour dessiner au hasard
np.random.seed(0)

n = 3 #Nombre de données de dessin
lower_indices = np.random.choice(np.arange(len(data_low)), n)
upper_indices = np.random.choice(np.arange(len(data_up)), n)
indices = np.hstack([lower_indices, upper_indices])

fig, axes = plt.subplots(ncols=n*2, nrows=2, figsize=(15, 15*2/(n*2)))

#Dessin de données avec une petite perte
for i in range(n):
    inputs = data_low[indices[i]][0]
    outputs = data_low[indices[i]][1]
    loss = data_low[indices[i]][2]
    axes[0, i].imshow(inputs.reshape(28, 28), cmap='gray')
    axes[1, i].imshow(outputs.reshape(28, 28), cmap='gray')
    axes[0, i].set_title('Input')
    axes[1, i].set_title('Small Loss={:.2f}'.format(loss))
    axes[0, i].set_xticks([])
    axes[0, i].set_yticks([])
    axes[1, i].set_xticks([])
    axes[1, i].set_yticks([])

#Dessin de données avec une grande perte
for i in range(n, 2*n):
    inputs = data_up[indices[i]][0]
    outputs = data_up[indices[i]][1]
    loss = data_up[indices[i]][2]
    axes[0, i].imshow(inputs.reshape(28, 28), cmap='gray')
    axes[1, i].imshow(outputs.reshape(28, 28), cmap='gray')
    axes[0, i].set_title('Input')
    axes[1, i].set_title('Large Loss={:.2f}'.format(loss))
    axes[0, i].set_xticks([])
    axes[0, i].set_yticks([])
    axes[1, i].set_xticks([])
    axes[1, i].set_yticks([])
plt.tight_layout()

Les images avec une perte faible ont tendance à avoir de grandes zones noires, tandis que les images avec une perte élevée ont tendance à avoir de grandes zones blanches. Puisque MNIST est fondamentalement une image avec une grande zone noire, si la zone blanche est grande, elle sera "étrangère" pour AE et il semble que la reconstruction ne réussisse pas.

Distance Approach k-NN(k-Nearest Neighbor) k-NN est une méthode pour mesurer le degré d'anomalie basé sur la distance entre les données. Dans k-NN, considérons un cercle contenant k points de données à proximité. Le chiffre ci-dessous est pour $ k = 5 $.

Comme le montre la figure, k-NN définit le rayon du cercle comme le degré d'anomalie basé sur l'idée que ** le cercle dessiné par les données anormales doit être plus grand que les données normales. Cependant, le seuil du rayon qui sépare normal / anormal du nombre de points voisins $ k $ doit être réglé de manière appropriée.

Faisons une expérience numérique. Le seuil est le point de division à 90% dans la distribution des valeurs de rayon des données normales.

python

#Créer des données
np.random.seed(123)

#Moyenne et variance
mean = np.array([2, 2])
cov = np.array([[3, 0], [0, 3]])

#Génération de données normales
norm = np.random.multivariate_normal(mean1, cov, size=50)

#Génération de données anormales(Généré à partir d'une distribution uniforme)
lower, upper = -10, 10
anom = (upper - lower)*np.random.rand(20,20) + lower

s=8
plt.scatter(norm[:,0], norm[:,1], s=s, color='green', label='normal')
plt.scatter(anom[:,0], anom[:,1], s=s, color='red', label='anomaly')
plt.legend()

Je ne savais pas comment terminer la détection d'anomalie par k-NN avec la bibliothèque, alors je l'ai fait moi-même. .. ..

python

#Une fonction qui calcule le rayon du cercle près de k de toutes les données
def kNN(data, k=5):
    """
    input: data:list, k:int
    output:Rayon des cercles près de k dans toutes les données:list
    
La distance adopte la distance euclidienne.
    """
    output = []
    for d in data:
        distances = []
        for p in data:
            distances.append(np.linalg.norm(d - p))
        output.append(sorted(distances)[k])
    return output

#1 si supérieur au seuil,Une fonction qui renvoie 0 si elle est petite
def kNN_scores(point, data, thres, K):
    distances = []
    for p in data:
        distances.append(np.linalg.norm(point - p))
    dist = sorted(distances)[K]
    if dist > thres:
        return 1
    else:
        return 0

python

def draw_KNN(norm, anom, thres, K, only_norm=False):
    # plot the line, the points, and the nearest vectors to the plane
    xx, yy = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
    Z = [kNN_scores(point, norm, thres, K) for point in np.c_[xx.ravel(), yy.ravel()]]
    Z = np.array(Z).reshape(xx.shape)

    plt.title("KNN:K="+str(K))
    plt.contourf(xx, yy, Z)
    
    if not only_norm:
        nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
        am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

        plt.axis('tight')
        plt.xlim((-10, 10))
        plt.ylim((-10, 10))
        plt.legend([nm, am],
                   ["normal", "anomaly"],
                   loc="upper left",
                   prop=matplotlib.font_manager.FontProperties(size=11))

    if only_norm:
        nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)

        plt.axis('tight')
        plt.xlim((-10, 10))
        plt.ylim((-10, 10))
        plt.legend([nm],
                   ["normal"],
                   loc="upper left",
                   prop=matplotlib.font_manager.FontProperties(size=11))
        
    plt.show()

python

for k in [1,5,30]:
    anomaly_scores = kNN(norm, k=k)
    thres = pd.Series(anomaly_scores).quantile(0.90)
    draw_KNN(norm, anom, thres, k)

Vous pouvez voir que lorsque k augmente, la limite normale / anormale devient plus lisse.

Maintenant, expérimentons avec un exemple légèrement désordonné. Supposons que ** les données normales soient générées à partir de deux distributions normales avec des densités différentes ** comme indiqué ci-dessous.

python

np.random.seed(123)

#Moyenne et variance
mean1 = np.array([3, 3])
mean2 = np.array([-5,-5])
cov1 = np.array([[0.5, 0], [0, 0.5]])
cov2 = np.array([[3, 0], [0, 3]])

#Génération de données normales(Généré à partir de deux distributions normales)
norm1 = np.random.multivariate_normal(mean1, cov1, size=50)
norm2 = np.random.multivariate_normal(mean2, cov2, size=10)
norm = np.vstack([norm1, norm2])

s=8
plt.scatter(norm[:,0], norm[:,1], s=s, color='green', label='normal')
plt.legend()

python

for k in [1,5,30]:
    anomaly_scores = kNN(norm, k=k)
    thres = pd.Series(anomaly_scores).quantile(0.90)
    draw_KNN(norm, [], thres, k, only_norm=True)

En vous concentrant particulièrement sur le cas de k = 1, vous pouvez voir que kNN ne fonctionne pas bien lorsque des données normales sont générées à partir de différents clusters avec des densités différentes. ** On peut imaginer qu'il peut être traité en ajoutant le processus «important l'écart par rapport à l'amas dense et tolérer l'écart par rapport à l'amas clairsemé» **.

En fait, c'est la méthode appelée LOF présentée dans la section suivante.

LOF (facteur de valeur aberrante locale)

Le cas où k-NN introduit dans la section précédente ne fonctionne pas bien est montré dans la figure.

Ce sera. Bien sûr, je voudrais augmenter le degré d'anomalie sur le côté gauche, mais LOF introduit le concept de ** densité de données ** et le réalise. Je vais omettre les détails mathématiques, mais je fais ce qui suit.

Le degré d'anomalie est défini par ** le rapport entre le degré de squishyness de soi-même et le degré de squishyness des points voisins **.

Expérimentons maintenant avec LOF.

python

def draw_LOF(norm, anom, K, only_norm=False):
    # plot the line, the points, and the nearest vectors to the plane
    xx, yy = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
    Z = - model.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = np.array(Z).reshape(xx.shape)

    plt.title("LOF:K="+str(K))
    plt.contourf(xx, yy, Z)
    
    if not only_norm:
        nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)
        am = plt.scatter(anom[:, 0], anom[:, 1], c='red', s=s)

        plt.axis('tight')
        plt.xlim((-10, 10))
        plt.ylim((-10, 10))
        plt.legend([nm, am],
                   ["normal", "anomaly"],
                   loc="upper left",
                   prop=matplotlib.font_manager.FontProperties(size=11))

    if only_norm:
        nm = plt.scatter(norm[:, 0], norm[:, 1], c='green', s=s)

        plt.axis('tight')
        plt.xlim((-10, 10))
        plt.ylim((-10, 10))
        plt.legend([nm],
                   ["normal"],
                   loc="upper left",
                   prop=matplotlib.font_manager.FontProperties(size=11))
        
    plt.show()

python

from sklearn.neighbors import LocalOutlierFactor

for k in [1, 5, 30]:
    model = LocalOutlierFactor(n_neighbors=k, novelty=True)
    model.fit(norm)
    draw_LOF(norm, [], k, only_norm=True)

De cette manière, il semble que la détection d'anomalies puisse être effectuée après avoir pris en compte la densité du cluster.

Problème de détection d'anomalies en évolution

Alors que j'étudie la détection d'anomalies à l'avenir, j'aimerais étudier les sujets suivants:

--Détection d'anomalies des données de séries chronologiques --Détection d'anomalies des données graphiques --Peut être utilisé pour la détection du blanchiment d'argent