Python离群值检测算法 -- Isolate Forest

什么是Isolate Forest？

许多离群点检测方法通常先分析正常数据点，然后找出不符合正常数据模式的观测值。然而，Liu、Ting和Zhou（2008）提出的Isolate Forest（IForest）与这些方法不同。相反，IForest直接识别异常点，而不是通过分析正常数据点来发现异常值。它使用树形结构来隔离每个观测点，异常点往往是最先被挑出来的数据点，而正常点则隐藏在树的深处。他们将每棵树称为Isolate Tree（iTree），构建了一个iTrees树群。异常点是指iTrees上平均路径长度较短的观测点。

iTree使用分区图和树来解释如何隔离数据点。红点最远离其他点，然后是绿点，最后是蓝点。在分区图中，只需一个 "切口 "就能将红点与其他点分开。第二个切点是绿点，第三个切点是蓝点，依此类推。分离一个点所需的切割次数越多，该点在树中的位置就越深。切割次数的倒数就是异常得分。树状结构也说明了同样的问题。通过一次分叉可以找到红点，第二次分叉可以找到绿点，第三次分叉可以找到蓝点，以此类推。深度数可以很好地代表异常点的得分。为了与异常点得分高的惯例保持一致，异常点得分被定义为深度数的倒数。

iTree是一种二叉树，每个节点都有0或2个子节点。树生长的条件包括：末端节点只有一个数据点，节点中的所有数据值相同，或者树达到研究人员设置的高度限制。iTree在所有末端节点都有一个数据点之前并不需要完全发展。通常情况下，当高度深度达到设定的限制时，树就会停止生长，因为我们关注的是靠近根节点的异常点。因此，构建一个大的iTree并不是必要的，因为iTree中的大部分数据都是正常数据点。小样本量能产生更好的iTree，因为沼泽效应和掩蔽效应会减弱。iTree算法与决策树算法不同，因为iTree不使用目标变量来训练树，它是一种无监督学习方法。

为什么是 "森林"

随机森林（Random Forests）可能比 "孤立森林"（Isolated Forests）更为常见。森林是指用于构建树木集合学习的概念。为何需要这样做呢？众所周知，单一决策树存在过拟合的缺点，这意味着模型对训练数据的预测效果很好，但对新数据的泛化效果较差。集合策略通过构建多棵决策树，然后对它们的预测结果进行平均，从而克服了这一问题。

由于孤立森林不使用任何距离度量来检测异常点，因此速度快，占用内存少。这一优势使其适用于大数据量和高维问题。

图 (B) 显示了一个数据矩阵，每一行都是一个具有多维值的观测值。IForest 的目标是为每个观测值分配离群值。首先，它会随机选择任意数量的行和任意数量的列来创建表格，如 (1)、(2) 和 (3)。一个观测值至少会出现在一个表格中。每个表格都会建立一棵 iTree 树，以显示离群点得分。表（1）有 6 行 3 列。表（1）的第一个切分点可能是第 6 个观测值，因为它的值与其他观测值非常不同。之后，表（1）的第二个切分点可能是第 4 个观测值。同样，在表（3）中，第一个切分点可能是第 6 个观测值（即第三条记录）。第二个切分点是第 4 个观测点（即表中的第一条记录）。简而言之，如果有N张表，就会有N个 iTrees。一个观测值最多可以有 N 个分数。IForest 会计算分数的算术平均值，得出最终分数。

建模流程

步骤 1：建立模型

我生成了一个包含六个变量和 500 个观测值的模拟数据集。无监督模型只使用 X 变量，而 Y 变量仅用于验证。异常值的百分比设定为 5%，contamination=0.05。可以绘制前两个变量的散点图，黄色的点表示异常值，紫色的点为正常数据点。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pyod.utils.data import generate_data
contamination = 0.05 # percentage of outliers
n_train = 500       # number of training points
n_test = 500        # number of testing points
n_features = 6      # number of features
X_train, X_test, y_train, y_test = generate_data(
    n_train=n_train, 
    n_test=n_test, 
    n_features= n_features, 
    contamination=contamination, 
    random_state=123)

X_train_pd = pd.DataFrame(X_train)
X_train_pd.head()

将树的大小 max_samples 设置为 40 个观测值。在 IForest 中，较小的样本量可以生成更好的 iTrees，无需指定较大的树规模。

from pyod.models.iforest import IForest
isft = IForest(contamination=0.05, max_samples=40, behaviour='new') 
isft.fit(X_train)

# Training data
y_train_scores = isft.decision_function(X_train)
y_train_pred = isft.predict(X_train)

# Test data
y_test_scores = isft.decision_function(X_test)
y_test_pred = isft.predict(X_test) # outlier labels (0 or 1)

# Threshold for the defined comtanimation rate
print("The threshold for the defined contamination rate:" , isft.threshold_)

def count_stat(vector):
    # Because it is '0' and '1', we can run a count statistic. 
    unique, counts = np.unique(vector, return_counts=True)
    return dict(zip(unique, counts))

print("The training data:", count_stat(y_train_pred))
print("The training data:", count_stat(y_test_pred))

The threshold for the defined contamination rate: 
-3.986394547794703e-15
The training data: {0: 475, 1: 25}
The training data: {0: 470, 1: 30}

污染率(contamination)

在实际应用中，通常我们无法确定异常值的百分比。在第 (C.2) 节中会说明，当我们事先无法确定异常值的百分比时，如何确定一个合理的阈值。PyOD 默认的污染率为 10%。在这里，我将污染率设置为 5%，因为在训练样本中污染率为 5%。这个参数不会影响离群值分数的计算。内置函数threshold_会根据污染率计算训练数据的阈值。在本例中，当污染率为 0.05 时，阈值为-5.082e-15。函数decision_functions()用来生成离群值，函数predict()则根据阈值分配标签（1 或 0）。

超参数

我将用.get_params() 来解释一些重要参数：

isft.get_params()

{'behaviour': 'new',
 'bootstrap': False,
 'contamination': 0.05,
 'max_features': 1.0,
 'max_samples': 40,
 'n_estimators': 100,
 'n_jobs': 1,
 'random_state': None,
 'verbose': 0}

• "max_samples"（最大样本数）指定训练每个基本估计子所提取的样本数量，这是一个重要参数；
• "n_estimators"表示集合中树的数量，默认为 100；
• "max_features"（最大特征）指定训练每个基本估计器所提取的特征数量，默认为 1.0；
• "n_jobs"表示并行运行的作业数量，如果设置为-1，则作业数将根据内核数确定。

特征重要性

IForest使用树形结构，能够衡量特征在确定异常值时的相对重要性，通过吉尼杂质指数来衡量特征的重要性，总和为1.0。

isft_vi = isft.feature_importances_
isft_vi

array([0.17325335, 0.13573292, 0.17200832,
       0.17157248, 0.17091259, 0.17652033])

可以像树模型一样绘制特征重要性图。下图显示了特征在确定异常值时的相对强度。

    
from matplotlib import pyplot as plt
for_plot = pd.DataFrame({'x_axis':X_train_pd.columns,
              'y_axis':isft_vi}).sort_values(by='y_axis',ascending=True)
for_plot['y_axis'].plot.barh()

步骤 2 - 确定模型的合理阈值

阈值应根据离群值的直方图来确定，下图建议阈值为0.0左右，这意味着大部分正常数据的离群值小于0.0，异常数据的离群值则处于较高范围。

import matplotlib.pyplot as plt
plt.hist(y_train_scores, bins='auto') # arguments are passed to np.histogram
plt.title("Outlier score")
plt.show()

第 3 步 - 显示正常组和异常组的描述性统计结果

正常组和异常组分析是验证模型合理性的关键步骤。我开发了一个名为descriptive_stat_threshold()的简要函数，用于展示正常组和异常组特征的大小和描述统计信息。此外，我还对污染率阈值进行了简单设置，您可以测试不同的阈值来确定合理的异常组大小。

threshold = isft.threshold_ # Or other value from the above histogram

def descriptive_stat_threshold(df,pred_score, threshold):
    # Let's see how many '0's and '1's.
    df = pd.DataFrame(df)
    df['Anomaly_Score'] = pred_score
    df['Group'] = np.where(df['Anomaly_Score']'Normal', 'Outlier')

    # Now let's show the summary statistics:
    cnt = df.groupby('Group')['Anomaly_Score'].count().reset_index().rename(columns={'Anomaly_Score':'Count'})
    cnt['Count %'] = (cnt['Count'] / cnt['Count'].sum()) * 100 # The count and count %
    stat = df.groupby('Group').mean().round(2).reset_index() # The avg.
    stat = cnt.merge(stat, left_on='Group',right_on='Group') # Put the count and the avg. together
    return (stat)

descriptive_stat_threshold(X_train,y_train_scores, threshold)

上表包含模型评估和结果的要点。特别注意要用特征名称标记特征，以有效展示。

• 离群组的大小： 离群值组的大小取决于所选的阈值。较高的阈值会使得该组规模较小。
• 每组中的特征统计数据： 特征统计数据应该与先前的业务知识一致。如果某些特征显示出令人费解的结果，应重新检查或删除该特征。需要反复模型以确保所有特征都是合理的。
• 平均异常值： 离群组的平均离群分数远高于正常组（0.18>-0.10）。不必对分数作过多解释。

由于我们已经拥有了基本事实而生成了数据，因此可以生成混淆矩阵来了解模型的性能。该模型表现出色，识别出了所有 25 个异常值。

def confusion_matrix(actual,score, threshold):
    Actual_pred = pd.DataFrame({'Actual': actual, 'Pred': score})
    Actual_pred['Pred'] = np.where(Actual_pred['Pred']<=threshold,0,1)
    cm = pd.crosstab(Actual_pred['Actual'],Actual_pred['Pred'])
    return (cm)
confusion_matrix(y_train,y_train_scores,threshold)

通过聚合多个模型实现模型稳定性

IForest算法对异常值非常敏感，可能会导致过拟合。为了得到稳定的预测结果，可以汇总多个模型的得分。在所有超参数中，树的数量n_estimators可能是最关键的参数。我会根据树的数量范围创建5个模型，然后取这些模型的平均预测值作为最终的模型预测值。PyOD模块提供了四种汇总结果的方法，您只需使用一种方法来生成汇总结果。请注意为这些函数安装pip install combo。

from pyod.models.combination import aom, moa, average, maximization
from pyod.utils.utility import standardizer
from pyod.models.iforest import IForest

# Standardize data
X_train_norm, X_test_norm = standardizer(X_train, X_test)

# Test a range of the number of trees
k_list = [100, 200, 300, 400, 500]
n_clf = len(k_list)
# Just prepare data frames so we can store the model results
train_scores = np.zeros([X_train.shape[0], n_clf])
test_scores = np.zeros([X_test.shape[0], n_clf])

# Modeling
for i in range(n_clf):
    k = k_list[i]
    #isft = IForest(contamination=0.05, max_samples=k) 
    isft = IForest(contamination=0.05, n_estimators=k) 
    isft.fit(X_train_norm)
    
    # Store the results in each column:
    train_scores[:, i] = isft.decision_function(X_train_norm) 
    test_scores[:, i] = isft.decision_function(X_test_norm) 
# Decision scores have to be normalized before combination
train_scores_norm, test_scores_norm = standardizer(train_scores,test_scores)

对 5 个模型的预测取平均值，得到离群值的平均值（"y_by_average"）。下图中绘制直方图。

# Combination by average
# The test_scores_norm is 500 x 10. The "average" function will take the average of the 10 columns. The result "y_by_average" is a single column: 
y_train_by_average = average(train_scores_norm)
y_test_by_average = average(test_scores_norm)
import matplotlib.pyplot as plt
plt.hist(y_train_by_average, bins='auto') # arguments are passed to np.histogram
plt.title("Combination by average")
plt.show()

上图表明阈值等于 1.0。因此，在下表中列出了正常组和离群组的特征。其中确定 25 个数据点为异常值。

descriptive_stat_threshold(X_train,
                           y_train_by_average,
                           1.0)

总结

大多数基于模型的异常检测方法都是构建正常实例的轮廓，然后将不符合正常轮廓的实例识别为异常值。相比之下，IForest 能直接、明确地隔离异常数据。IForest 采用树形结构来隔离每一个数据点，异常点被首先挑出，而正常点则往往聚集在树状结构中。由于Isolation Forest 不使用任何距离度量来检测异常点，因此速度快，适用于大数据量和高维问题。