Python 之 异常检测/异常点搜索

 

异常点搜索

继续话题：对离群值的处理。

已经开发了一套Python工具，用于检测数据中的排放和搜索异常点--Python异常点检测(PyOD)

PyOD是一个全面的、可扩展的Python工具包，用于检测多维数据中的外部对象（排放）。

PyOD的优势。

1. 统一的API，详细的文档和不同算法的交互式例子。
2. 高级模型，包括神经网络/深度训练和发射集合。
3.使用JIT优化性能，并在可能的情况下使用numba进行并行。
和joblib。
4. 同时兼容Python 2和Python 3。
协议：线性模型、基于近似的、概率算法、离群组、神经网络。

算法

1.导入所需模块

import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

import numpy as np
import pandas as pd
from pyod.utils.data import generate_data, get_outliers_inliers

from scipy import stats

import itertools

itertools模块实现了创建迭代器以有效组织循环的函数。

2.使用函数generate_data创建排放量为10%的数据

contamination = 0.1 # percentage of outliers
n_train = 500       # number of training points
n_test = 500        # number of testing points
n_features = 2      # number of features
X_train, y_train, X_test, y_test = generate_data(n_train=n_train, 
                                                 n_test=n_test, 
                                                 n_features=n_features, 
                                                 contamination=contamination)

你可以分别计算离群值（不排放）和离群值（排放）的数量。

# store outliers and inliers in different numpy arrays
x_outliers, x_inliers = get_outliers_inliers(X_train, y_train)

n_inliers = len(x_inliers)
n_outliers = len(x_outliers)

print('inliers: ', n_inliers, ' outliers: ', n_outliers)

3.在散点图上显示数据。

# Plot
plt.scatter(X_train[:, 0], X_train[:, 1], alpha=0.8)
plt.title('Scatter plot')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

图中清晰地显示了云--一组离群者（不是排放）和云附近的少量排放点。

我们使用算法k最近的邻居（k-NN）来识别异常。

k-NN算法是一种非参数化的方法，可以识别出最接近的教程例子𝑘。任何孤立的数据点都有可能被归类为排放。

以下代码行为k-NN模型提供训练，并将模型存储为clf。

# train kNN detector
from pyod.models.knn import KNN
clf_name = 'KNN'
clf = KNN()
clf.fit(X_train)
# If you want to see the predictions of the training data, you can use this way:
y_train_scores = clf.decision_scores_

让我们在图上显示教程数据集的收到的 "分数"（估计每个数据点有多少是离群值或不是）。

plt.plot(y_train_scores, linestyle='--', color='blue', linewidth=3)
plt.plot(y_train, linestyle='-', color='red', linewidth=3)
plt.legend(['estimated scores on train data', 'set for train data'])
plt.show();

通过训练好的k-NN模型，你可以预测一个点是否是测试数据集的离群值。decision_functions函数根据模型预测生成对异常情况的估计。

# Now we have the trained K-NN model, let's apply to the test data to get the predictions
y_test_pred = clf.predict(X_test) # outlier labels (0 or 1)

让我们绘制数据集的获得的估计值。

plt.plot(y_test_pred, linestyle='--', color='blue', linewidth=3)
plt.plot(y_test, linestyle='-', color='red', linewidth=3)
plt.legend(['estimated value on test data', 'set for test data'])
plt.show();

我们将计算模型所估计的排放量百分比。

 

# 因为是'0'和'1'，我们可以进行一个计数统计。异常的数量大概是百分之十，我们之前已经产生了。
unique, counts = np.unique(y_test_pred, return_counts=True)
result = dict(zip(unique, counts))
print(result)
print(100 * result[1] / (result[0] + result[1]), '% of anomalies')

模型估计的排放百分比为10%（我们原来设定的排放百分比）。

让我们绘制测试样本的 "分数"（对 "异常 "观测值的估计）。

# And you can generate the anomaly score using clf.decision_function:
y_test_scores = clf.decision_function(X_test)

plt.plot(y_test_scores, linestyle='--', color='blue', linewidth=3)
plt.plot(y_test, linestyle='-', color='red', linewidth=3)
plt.legend(['estimated scores on test data', 'set for test data'])
plt.show();

你可以使用mlxtend包将k-NN算法的结果可视化。

Mlxtend（机器学习扩展）是一个Python工具库，用于日常数据处理任务。

from mlxtend.plotting import plot_decision_regions

# Plotting decision region
plot_decision_regions(X_train, y_train.astype('int'), clf=clf, legend=2)  # y only as int variables
# Adding axes annotations
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Knn decision regions results')
plt.show();

如何界定 "异常评估"？回想一下，k-NN模型使用欧几里得距离来衡量距离。排放是相邻点的远距离点，所以排放估计是由距离值决定的。我们的任务是找到高离群值的点。我们可以用直方图来寻找这些点。

plt.hist(y_test_scores, bins='auto')  # arguments are passed to np.histogram
plt.title("Histogram with 'auto' bins")
plt.show();

此外，您还可以计算正常和异常聚类的统计。

在我们的案例中，我们有两个特征（标记为 "0 "和 "1"）。第1群组的平均异常估计值远高于第0群组，汇总统计也显示出这两个群组之间的显著差异。因此，我们认为群组1中的数据点可能存在异常，值得进一步分析。

# Let's see how many '0's and '1's.
df_test = pd.DataFrame(X_test)
df_test['score'] = y_test_scores
df_test['cluster'] = np.where(df_test['score'] < 1, 0, 1)
df_test['cluster'].value_counts()

# Now let's show the summary statistics:
df_test.groupby('score').mean()
df_test.groupby('cluster').mean()

4.在pyod的帮助下，可以同时实现多个模型来搜索排放。

让我们在分类器字典中设置所选模型。

from pyod.models.abod import ABOD

classifiers = {
     'Angle-based Outlier Detector (ABOD)' : ABOD(contamination=contamination),
     'K Nearest Neighbors (KNN)' :  KNN(contamination=contamination)
}

让我们创建一个网格显示可视化结果。

# create a meshgrid 
xx, yy = np.meshgrid(np.linspace(-10, 10, 200), np.linspace(-10, 10, 200))

以下代码行是按顺序保护所选模型，以及展示可视化所获得的结果。

# set the figure size
plt.figure(figsize=(10, 10))

for i, (clf_name, clf) in enumerate(classifiers.items()) :
    # fit the dataset to the model
    clf.fit(X_train)

    # predict raw anomaly score
    scores_pred = clf.decision_function(X_train)*-1

    # prediction of a datapoint category outlier or inlier
    y_pred = clf.predict(X_train)

    # no of errors in prediction
    n_errors = (y_pred != y_train).sum()
    print('No of Errors : ', clf_name, n_errors)

    # rest of the code is to create the visualization

    # threshold value to consider a datapoint inlier or outlier
    threshold = stats.scoreatpercentile(scores_pred, 100 * contamination)

    # decision function calculates the raw anomaly score for every point
    Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()]) * -1
    Z = Z.reshape(xx.shape)

    subplot = plt.subplot(1, 2, i + 1)

    # fill blue colormap from minimum anomaly score to threshold value
    subplot.contourf(xx, yy, Z, levels = np.linspace(Z.min(), threshold, 10),cmap=plt.cm.Blues_r)

    # draw red contour line where anomaly score is equal to threshold
    a = subplot.contour(xx, yy, Z, levels=[threshold], linewidths=2, colors='red')

    # fill orange contour lines where range of anomaly score is from threshold to maximum anomaly score
    subplot.contourf(xx, yy, Z, levels=[threshold, Z.max()], colors='orange')

    # scatter plot of inliers with white dots
    b = subplot.scatter(X_train[:-n_outliers, 0], X_train[:-n_outliers, 1], c='white', s=20, edgecolor='k') 
    # scatter plot of outliers with black dots
    c = subplot.scatter(X_train[-n_outliers:, 0], X_train[-n_outliers:, 1], c='black', s=20, edgecolor='k')
    subplot.axis('tight')

    subplot.legend(
        [a.collections[0], b, c],
        ['learned decision function', 'true inliers', 'true outliers'],
        prop=matplotlib.font_manager.FontProperties(size=10),
        loc='lower right')

    subplot.set_title(clf_name)
    subplot.set_xlim((-10, 10))
    subplot.set_ylim((-10, 10))
    
plt.show(); 

5.实验

在有些情况下，需要在一行中找到一个点，从这个点开始，所研究的过程会急剧改变其特征。

例如，在工业厂房中，需要根据测量仪器的读数提前（尽可能多的时间）预测紧急情况。

data_with_noise = pd.read_csv('data/samples/add_noise_1_7603.csv')
data_with_trend = pd.read_csv('data/samples/change_trend_10_7945.csv')
data_with_decreased_dispersion = pd.read_csv('data/samples/decrease_dispersion_190_8980.csv')
data_with_increased_dispersion = pd.read_csv('data/samples/increase_dispersion_190_9497.csv')
data_with_shifted_trend = pd.read_csv('data/samples/shift_trend_10_8709.csv')

data = pd.concat([data_with_noise, data_with_trend, data_with_shifted_trend, 
                  data_with_decreased_dispersion, data_with_increased_dispersion], axis=1)
data.columns = ['with_noise', 'change_trend', 'shift_trend', 'decrease_dispersion', 'increase_dispersion']
points = [7603, 7945, 8980, 9497, 8709]

智能化数据

data['with_noise'].plot();

增加趋势数据

data['change_trend'].plot();

变化趋势数据

data['shift_trend'].plot();

降低数据的分散度

data['decrease_dispersion'].plot();

增加数据的分散性

data['increase_dispersion'].plot();

让我们试着确保k-NN。

fig, axes = plt.subplots(nrows=data.shape[1], ncols=2, figsize=(10,10))

for idx, column in enumerate(data.columns):
    print(column, ' is processed')
    print('anomaly point is ', points[idx])
    clf_tsu_data = KNN()
    clf_tsu_data.fit(data[[column]])
    scores = clf_tsu_data.decision_scores_
    predictions = clf_tsu_data.predict(data[[column]])
    
    try:
        print('found anomaly point is ', list(predictions[int(data.shape[0] / 2):]).index(1) + int(data.shape[0] / 2))
    except:
        print('no anomaly point')
    
    axes[idx, 0].plot(scores, linestyle='--', color='blue', linewidth=3)
    axes[idx, 0].set_title(column + '_scores')
    
    axes[idx, 1].plot(predictions, linestyle='-', color='red', linewidth=3)
    axes[idx, 1].set_title(column + '_predictions')
    
    fig.tight_layout();

 

让我们来考虑各种数据转换和其中的结构搜索吧

我们使用另一个Python库的工具来检测异常--luminol。

import luminol
from luminol.anomaly_detector import AnomalyDetector
from luminol.correlator import Correlator

6.取第一差值和滚动中位数

帮助寻找有噪声和变化趋势的数据的异常点。

diff_data = data.diff()
diff_data.dropna(inplace=True)

window = 100
diff_data = diff_data.rolling(window).median()
diff_data.dropna(inplace=True)

fig, axes = plt.subplots(nrows=data.shape[1], ncols=1, figsize=(5,10))
for idx, column in enumerate(diff_data.columns):
    axes[idx].plot(diff_data[column])
    axes[idx].set_title(column + '_diff_rolling_median')
    
    fig.tight_layout(); 

diff_data_dict = diff_data.to_dict()  

fig, axes = plt.subplots(nrows=data.shape[1], ncols=1, figsize=(5,10))

for idx, column in enumerate(diff_data.columns):
    detector = AnomalyDetector(diff_data_dict[column])
    anomalies = detector.get_anomalies()
    scores = detector.get_all_scores()
    
    axes[idx].plot(scores.values)
    axes[idx].set_title(column + '_scores')
    
    print(column, ' is processed')
    print('true anomaly point is ', points[idx])
    try:
        print('found anomaly point is ', next(x[0] for x in enumerate(scores.values) if x[1] > 0)) 
    except:
        print('no anomaly point')
    
    fig.tight_layout(); 

7.取第一个差值并滚动变异 rolling var

有助于找到噪声数据的异常点，提高分散度。

diff_data_var = data.diff()
diff_data_var.dropna(inplace=True)

window = 100
diff_data_var = diff_data_var.rolling(window).var()
diff_data_var.dropna(inplace=True)

fig, axes = plt.subplots(nrows=data.shape[1], ncols=1, figsize=(5,10))
for idx, column in enumerate(diff_data_var.columns):
    axes[idx].plot(diff_data_var[column])
    axes[idx].set_title(column + '_diff_rolling_var')
    
    fig.tight_layout(); 

 

diff_data_dict = diff_data_var.to_dict()  

fig, axes = plt.subplots(nrows=data.shape[1], ncols=1, figsize=(5,10))

for idx, column in enumerate(diff_data_var.columns):
    detector = AnomalyDetector(diff_data_dict[column])
    anomalies = detector.get_anomalies()
    scores = detector.get_all_scores()
    
    axes[idx].plot(scores.values)
    axes[idx].set_title(column + '_scores')
    
    print(column, ' is processed')
    print('true anomaly point is ', points[idx])
    try:
        print('found anomaly point is ', next(x[0] for x in enumerate(scores.values) if x[1] > 0)) 
    except:
        print('no anomaly point')
    
    fig.tight_layout(); 

8.取二次差值和滚动中位数

有助于发现分散度和移动趋势增加的数据的异常情况。

diff_data_ = data.diff()
diff_data_.dropna(inplace=True)

diff_data_ = diff_data_.diff()
diff_data_.dropna(inplace=True)

window = 100
diff_data_ = diff_data_var.rolling(window).var()
diff_data_.dropna(inplace=True)

ig, axes = plt.subplots(nrows=data.shape[1], ncols=1, figsize=(5,10))
for idx, column in enumerate(diff_data_.columns):
    axes[idx].plot(diff_data_[column])
    axes[idx].set_title(column + '_diff_diff')
    
    fig.tight_layout(); 

diff_data_dict = diff_data_.to_dict()  

fig, axes = plt.subplots(nrows=data.shape[1], ncols=1, figsize=(5,10))

for idx, column in enumerate(diff_data_.columns):
    detector = AnomalyDetector(diff_data_dict[column])
    anomalies = detector.get_anomalies()
    scores = detector.get_all_scores()
    
    axes[idx].plot(scores.values)
    axes[idx].set_title(column + '_scores')
    
    print(column, ' is processed')
    print('true anomaly point is ', points[idx])
    try:
        print('found anomaly point is ', next(x[0] for x in enumerate(scores.values) if x[1] > 0)) 
    except:
        print('no anomaly point')
    
    fig.tight_layout();