异常值检测方法比较——基于美国职业棒球联盟2023赛季击球数据

译者 | 朱先忠

审校 | 重楼

异常值检测是一种无监督的机器学习任务，用于识别给定数据集中的异常（即“异常观测”）。在大量现实世界中，当我们的可用数据集已经被异常“污染”时，异常值检测任务对于整个机器学习环节来说是非常有帮助的。当前，开源框架Scikit learn已经实现了几种异常值检测算法，在我们持有未受污染的基准数据的情况下，我们也可以使用这些算法进行新颖检测。这是一项半监督任务，可以预测新的观测结果是否为异常值。

概述

本文中，我们将比较的四种异常值检测算法如下：

• 椭圆包络（Elliptic Envelope）适用于低维正态分布数据。顾名思义，它使用多元正态分布来创建一个距离度量指标，以便将异常值与内部值区分开来。
• 局部异常因子（Local Outlier Factor）是一种针对观测的局部密度与其邻居的局部密度的比较。密度比其邻居低得多的观测被视为异常值。
• 具有随机梯度下降的单类支持向量机（One-Class Support Vector Machine with Stochastic Gradient Descent）是一类SVM的O（n）近似解。请注意，O（n²）数量级的单类SVM在我们的小示例数据集上运行良好，但对于您的实际应用场景来说可能不切实际。
• 孤立森林（Isolation Forest）是一种基于树的方法，其中异常值通过随机拆分方法比使用内部值方法可实现更快的隔离。

由于我们的任务是无监督的；所以，我们没有基本事实来比较这些算法的准确性。相反，我们想看看他们的结果（尤其是球员排名方面）彼此之间的差异，并对他们的行为和局限性获得一些直觉，这样我们就可以知道什么时候更喜欢其中一名球员而不是另一名。

接下来，让我们使用2023年美国职业棒球大联盟（MLB）赛季击球手表现的两个指标来比较其中的一些技术吧：

• 上垒率（OBP），击球手每次上垒时（通过击球、保送或投球命中）上垒的速度
• Slugging（SLG），每次击球的平均总垒数

击球手表现还有许多更复杂的指标，例如OBP加SLG（OPS）、基本平均加权（wOBA）和调整后的加权跑动（WRC+）。然而，我们将看到，除了常用和易于理解之外，OBP和SLG具有适度的相关性和近似正态分布，使它们非常适合这种比较。

数据集准备

我们使用pybalball包来获取击球数据。该Python包接受MIT许可，能够取得来自于Fangraphs.com、Baseball-Reference.com和其他来源的数据；当然，这些来源反过来又是从美国职业棒球大联盟获得的官方记录数据。

我们使用pybalball的2023年击球统计数据，可以通过batting_stats（FanGraphs网站）或batting_stasts_bref（权威网站Baseball Reference）获得。事实证明，FanGraphs中的球员名称格式更正确，但Baseball Reference中的球员团队和联盟在交易球员的情况下格式更好。为了获得一个具有更高可读性的数据集，我们实际上需要合并三个表：FanGraphs、Baseball Reference和键查询。

from pybaseball import (cache, batting_stats_bref, batting_stats, playerid_reverse_lookup) import pandas as pd cache.enable() # 重新运行时避免不必要的请求 MIN_PLATE_APPEARANCES = 200 # 为了可读性和合理的默认排序顺序 df_bref = batting_stats_bref(2023).query(f"PA >= {MIN_PLATE_APPEARANCES}" ).rename(columns={"Lev":"League", "Tm":"Team"} ) df_bref["League"] = \ df_bref["League"].str.replace("Maj-","").replace("AL,NL","NL/AL" ).astype('category') df_fg = batting_stats(2023, qual=MIN_PLATE_APPEARANCES) key_mapping = \ playerid_reverse_lookup(df_bref["mlbID"].to_list(), key_type='mlbam' )[["key_mlbam","key_fangraphs"] ].rename(columns={"key_mlbam":"mlbID", "key_fangraphs":"IDfg"} ) df = df_fg.drop(columns="Team" ).merge(key_mapping, how="inner", on="IDfg" ).merge(df_bref[["mlbID","League","Team"]], how="inner", on="mlbID" ).sort_values(["League","Team","Name"])

数据探索

首先，我们注意到这些指标的均值和方差是不同的，并且它们之间具有适度的相关性。我们还注意到，每个指标都是相当对称的，中值接近平均值。

print(df[["OBP","SLG"]].describe().round(3)) print(f"\nCorrelation: {df[['OBP','SLG']].corr()['SLG']['OBP']:.3f}")

让我们使用以下方法来可视化上面的联合分布：

• 球员散点图，全国职业棒球联盟（NL）对美国联盟（AL）着色
• 球员的二元核密度估计器（KDE）图，用高斯核平滑散射图以估计密度
• 每个指标的边际KDE图

import matplotlib.pyplot as plt import seaborn as sns g = sns.JointGrid(data=df, x="OBP", y="SLG", height=5) g = g.plot_joint(func=sns.scatterplot, data=df, hue="League", palette={"AL":"blue","NL":"maroon","NL/AL":"green"}, alpha=0.6 ) g.fig.suptitle("On-base percentage vs. Slugging\n2023 season, min " f"{MIN_PLATE_APPEARANCES} plate appearances" ) g.figure.subplots_adjust(top=0.9) sns.kdeplot(x=df["OBP"], color="orange", ax=g.ax_marg_x, alpha=0.5) sns.kdeplot(y=df["SLG"], color="orange", ax=g.ax_marg_y, alpha=0.5) sns.kdeplot(data=df, x="OBP", y="SLG", ax=g.ax_joint, color="orange", alpha=0.5 ) df_extremes = df[ df["OBP"].isin([df["OBP"].min(),df["OBP"].max()]) | df["OPS"].isin([df["OPS"].min(),df["OPS"].max()]) ] for _,row in df_extremes.iterrows(): g.ax_joint.annotate(row["Name"], (row["OBP"], row["SLG"]),size=6, xycoords='data', xytext=(-3, 0), textcoords='offset points', ha="right", alpha=0.7) plt.show()

散点图的右上角显示了一组优秀的击球，对应于SLG和OBP分布的上部重尾。这个小团体擅长上垒和多垒击球。我们认为他们在多大程度上是异常值（因为他们与大多数球员群体的距离），而不是内部值（因为他们彼此接近），这取决于我们选择的算法所使用的定义。

应用异常值检测算法

Scikit-learn的异常值检测算法通常提供了fit()和predict()两种方法，但算法之间也有例外，也有不同之处。我们将单独考虑每个算法，但我们将每个算法拟合到每个球员（m=453）的属性矩阵（n=2）中。然后，我们不仅会为每个球员打分，还会为每个属性范围内的值网格打分，以帮助我们可视化预测函数。

为了可视化决策边界，我们需要采取以下步骤：

1. 创建输入特征值的二维网格。
2. 将decision_function应用于网格上的每个点，这需要拆开网格。
3. 将预测重新成形为一个网格。
4. 绘制预测。

在此，我们将使用200x200大小的网格来覆盖现有的观测结果和一些填充，但您可以将网格调整到所需的速度和分辨率。

import numpy as np X = df[["OBP","SLG"]].to_numpy() GRID_RESOLUTION = 200 disp_x_range, disp_y_range = ( (.6*X[:,i].min(), 1.2*X[:,i].max()) for i in [0,1] ) xx, yy = np.meshgrid(np.linspace(*disp_x_range, GRID_RESOLUTION), np.linspace(*disp_y_range, GRID_RESOLUTION) ) grid_shape = xx.shape grid_unstacked = np.c_[xx.ravel(), yy.ravel()]

椭圆包络

椭圆包络的形状由数据的协方差矩阵确定，该协方差矩阵给出了特征i在主对角线[i，i]上的方差以及特征i和j在[i，j]位置上的协方差。由于协方差矩阵对异常值敏感，因此该算法使用最小协方差行列式（MCD）估计器，该估计器被推荐用于单峰和对称分布，通过随机混洗后的random_state输入参数以获得再现性。请注意，这个稳健的协方差矩阵稍后将再次派上用场。

因为我们想比较他们排名中的异常值分数，而不是二元异常值/内部分类，所以我们决定使用decision_function函数为球员打分。

from sklearn.covariance import EllipticEnvelope ell = EllipticEnvelope(random_state=17).fit(X) df["outlier_score_ell"] = ell.decision_function(X) Z_ell = ell.decision_function(grid_unstacked).reshape(grid_shape)

局部异常值因子

这种测量隔离度的方法基于k近邻（KNN）算法。我们计算每个观测到其最近邻居的总距离，以定义局部密度，然后将每个观测的局部密度与其邻近的局部密度进行比较。局部密度远小于其邻近的观测结果即被视为异常值。

选择要包括的邻近数量

K = int(np.sqrt(X.shape[0])) print(f"Using K={K} nearest neighbors.")

选择一个距离指标

在计算马氏距离时，我们将使用鲁棒性强的协方差矩阵。如果我们还没有通过椭圆包络计算它，那么我们可以直接计算它。

from scipy.spatial.distance import pdist, squareform # 如果我们还没有椭圆包络，我们可以计算鲁棒性强的协方差： # from sklearn.covariance import MinCovDet # robust_cov = MinCovDet().fit(X).covariance_ # 但我们可以从椭圆包络中重复使用它： robust_cov = ell.covariance_ print(f"Robust covariance matrix:\n{np.round(robust_cov,5)}\n") inv_robust_cov = np.linalg.inv(robust_cov) D_mahal = squareform(pdist(X, 'mahalanobis', VI=inv_robust_cov)) print(f"Mahalanobis distance matrix of size {D_mahal.shape}, " f"e.g.:\n{np.round(D_mahal[:5,:5],3)}...\n...\n")

拟合局部异常因子

还要注意，与其他算法不同，使用LOF时，我们被指示不要使用score_samples方法对现有观测值进行评分；该方法应仅用于新颖性检测。

from sklearn.neighbors import LocalOutlierFactor lof = LocalOutlierFactor(n_neighbors=K, metric="precomputed", novelty=True ).fit(D_mahal) df["outlier_score_lof"] = lof.negative_outlier_factor_

创建决策边界

from scipy.spatial.distance import cdist D_mahal_grid = cdist(XA=grid_unstacked, XB=X, metric='mahalanobis', VI=inv_robust_cov ) Z_lof = lof.decision_function(D_mahal_grid).reshape(grid_shape)

支持向量机（SGD-One-Class SVM）

SVM使用核技巧将特征转换到更高的维度，这样可以方便识别分离超平面。径向基函数（RBF）内核要求输入标准化，但正如StandardScaler的文档所指出的，该缩放器对异常值很敏感，所以我们将使用RobustScaler。正如SGDOneClassSVM的文档所建议的那样，我们将缩放后的输入管道传输到Nyström核近似中。

from sklearn.pipeline import make_pipeline from sklearn.preprocessing import RobustScaler from sklearn.kernel_approximation import Nystroem from sklearn.linear_model import SGDOneClassSVM suv = make_pipeline( RobustScaler(), Nystroem(random_state=17), SGDOneClassSVM(random_state=17) ).fit(X) df["outlier_score_svm"] = suv.decision_function(X) Z_svm = suv.decision_function(grid_unstacked).reshape(grid_shape)

孤立森林

这种基于树的隔离测量方法执行随机递归划分。如果隔离给定观测所需的平均分裂次数较低，则该观测被视为更强的候选异常值。与随机森林和其他基于树的模型一样，孤立森林并不假设特征是正态分布的，也不要求对其进行缩放。默认情况下，它构建100棵树。我们的示例仅使用两个特征，因此并没有启用特征采样。

from sklearn.ensemble import IsolationForest iso = IsolationForest(random_state=17).fit(X) df["outlier_score_iso"] = iso.score_samples(X) Z_iso = iso.decision_function(grid_unstacked).reshape(grid_shape)

结果：检查决策边界

请注意，来自上述这些模型的预测生成了不同的分布结果。我们选择应用QuantileTransformer转换器，目的是使预测值在给定网格上更具视觉可比性。从官方文档中，请注意如下提示：

这种变换是非线性的。它可能会扭曲在同一测量标准下测量的变量之间的线性相关性，但会使在不同测量标准上测量的变量更直接地具有可比性。

from adjustText import adjust_text from sklearn.preprocessing import QuantileTransformer N_QUANTILES = 8 # 每个图表有这么多颜色断点 N_CALLOUTS=15 # 为每个图表标记这么多顶部异常值 fig, axs = plt.subplots(2, 2, figsize=(12, 12), sharex=True, sharey=True) fig.suptitle("Comparison of Outlier Identification Algorithms",size=20) fig.supxlabel("On-Base Percentage (OBP)") fig.supylabel("Slugging (SLG)") ax_ell = axs[0,0] ax_lof = axs[0,1] ax_svm = axs[1,0] ax_iso = axs[1,1] model_abbrs = ["ell","iso","lof","svm"] qt = QuantileTransformer(n_quantiles=N_QUANTILES) for ax, nm, abbr, zz in zip( [ax_ell,ax_iso,ax_lof,ax_svm], ["Elliptic Envelope","Isolation Forest", "Local Outlier Factor","One-class SVM"], model_abbrs, [Z_ell,Z_iso,Z_lof,Z_svm] ): ax.title.set_text(nm) outlier_score_var_nm = f"outlier_score_{abbr}" qt.fit(np.sort(zz.reshape(-1,1))) zz_qtl = qt.transform(zz.reshape(-1,1)).reshape(zz.shape) cs = ax.contourf(xx, yy, zz_qtl, cmap=plt.cm.OrRd.reversed(), levels=np.linspace(0,1,N_QUANTILES) ) ax.scatter(X[:, 0], X[:, 1], s=20, c="b", edgecolor="k", alpha=0.5) df_callouts = df.sort_values(outlier_score_var_nm).head(N_CALLOUTS) texts = [ ax.text(row["OBP"], row["SLG"], row["Name"], c="b", size=9, alpha=1.0) for _,row in df_callouts.iterrows() ] adjust_text(texts, df_callouts["OBP"].values, df_callouts["SLG"].values, arrowprops=dict(arrowstyle='->', color="b", alpha=0.6), ax=ax ) plt.tight_layout(pad=2) plt.show() for var in ["OBP","SLG"]: df[f"Pctl_{var}"] = 100*(df[var].rank()/df[var].size).round(3) model_score_vars = [f"outlier_score_{nm}" for nm in model_abbrs] model_rank_vars = [f"Rank_{nm.upper()}" for nm in model_abbrs] df[model_rank_vars] = df[model_score_vars].rank(axis=0).astype(int) # Averaging the ranks is arbitrary; we just need a countdown order df["Rank_avg"] = df[model_rank_vars].mean(axis=1) print("Counting down to the greatest outlier...\n") print( df.sort_values("Rank_avg",ascending=False ).tail(N_CALLOUTS)[["Name","AB","PA","H","2B","3B", "HR","BB","HBP","SO","OBP", "Pctl_OBP","SLG","Pctl_SLG" ] + [f"Rank_{nm.upper()}" for nm in model_abbrs] ].to_string(index=False) )

分析与结论

看起来上面四种实现方案在如何定义异常值方面基本保持一致，但在分值和易用性方面存在一些明显差异。

椭圆包络方法

椭圆包络在没有配置的情况下也很容易使用，并且它提供了具有鲁棒性的协方差矩阵。如果您有低维数据，并且合理地期望它是正态分布的（通常情况并非如此），那么您可能想先尝试这种简单的方法。

单分类算法

值得注意的是，RBF内核需要一个额外的步骤来进行标准化，但在这个简单的例子中，它似乎也能很好地工作，而无需进行微调。

局部异常因子

另外，值得注意的是，LOF实现是最耗时的，因为我们要创建的是用于拟合和评分的稳健的距离矩阵。当然，我们本可以花一些时间来调整K。

孤立森林

与常见的基于监督树的学习器一样，隔离森林不进行外推分析，这使得它更适合于拟合受污染的数据集进行异常值检测，而不是拟合无异常的数据集用于新颖性检测（在这种情况下，它不会比现有观察更有力地对新的异常值进行评分）。

尽管“隔离森林”开箱即用效果很好，但它未能将大谷昭平列为棒球（可能是所有职业运动）中最大的异常值，这说明了任何异常值检测器的主要局限性：用于拟合它的数据。

需要说明的是，在我们的实验数据中，我们不仅省略了防守数据（对不起，奥斯汀·赫奇斯），也没有考虑投球数据。因为投手们甚至不再尝试击球了……除了Ohtani，他的赛季包括棒球史上第二好的打击率（BAA）和第11好的平均得分（ERA）（至少投了100局），一场完整的比赛被淘汰，以及一场他三振十名击球手并击出两支全垒打的比赛。

有人认为大谷昭平是一个模仿人类的高级外星人，但似乎更有可能有两个模仿同一个人的高级外星人。不幸的是，其中一人刚刚做了肘部手术，2024年不会投球了……但另一人刚刚签署了一份创纪录的10年7亿美元的合同。哈哈，得益于异常值检测，现在我们终于可以明白其中有关的原因了！

译者介绍

朱先忠，51CTO社区编辑，51CTO专家博客、讲师，潍坊一所高校计算机教师，自由编程界老兵一枚。

原文标题：