比较离群值检测方法
离群值检测方法比较
使用2023年美国职业棒球大联盟的击球数据
异常检测是一项无监督的机器学习任务，用于识别给定数据集中的异常（异常观察）。在许多现实世界的情况下，我们的可用数据集已经被异常“污染”。Scikit-learn 实现了几种异常检测算法，在我们有一个不受异常干扰的基线的情况下，我们也可以使用这些算法进行新颖性检测，这是一种半监督任务，可以预测新观测是否为异常值。
概述
我们将比较四种异常检测算法：
- 椭圆包络适用于低维正态分布的数据。顾名思义，它使用多元正态分布来创建一个距离度量，以将异常值与非异常值分开。
- 局部异常因子是一种将观察值的局部密度与其邻居的密度进行比较的方法。与其邻居相比，密度明显较低的观察值被视为异常值。
- 带有随机梯度下降（SGD）的单类支持向量机（SVM）是单类SVM的O(n)近似解。请注意，O(n²)的单类SVM在我们的小范例数据集上工作良好，但在实际使用中可能不切实际。
- 孤立森林是一种基于树的方法，其中异常值比非异常值更快地被随机分割孤立。
由于我们的任务是无监督的，我们没有基准来比较这些算法的准确性。相反，我们希望看到它们的结果（尤其是球员排名）在彼此之间的差异，并对其行为和限制有一些直观的认识，以便我们可以知道何时更喜欢其中一个。
让我们使用2023年美国职业棒球大联盟（MLB）赛季的击球数据来比较这些技术的几个指标：
- 大型语言模型（LLMs）的软硬件协同优化策略（SW/HW Co-optimization Strategy）
- 如何使用Langchain？逐步指南
- 这篇 AI 论文揭示了 X-Raydar：开创性的开源深度神经网络用于胸部 X 射线异常检测
- 上垒率（OBP），即击球手每次上场击球（击球、保送或被击球打中）时达到垒的率
- 嘲笑（SLG），即每场打席平均的垒打数
还有许多更复杂的击球表现指标，包括OBP加SLG（OPS），调整的上垒平均（wOBA）和调整权重得分（WRC+）。然而，我们将看到，除了被广泛使用和易于理解之外，OBP和SLG也具有适度的相关性和近似正态分布，使它们非常适合进行比较。
数据集准备
我们使用pybaseball包获取击球数据。这个Python包根据MIT 许可证返回来自Fangraphs.com、Baseball-Reference.com和其他来源的数据，这些来源又从美国职业棒球大联盟获得了官方记录。
我们使用pybaseball的2023年击球统计数据，可以通过batting_stats（FanGraphs）或batting_stats_bref（Baseball Reference）获取。事实证明，从Fangraphs获取的球员姓名格式更正确，但是在交易球员的情况下，来自Baseball Reference的球员队伍和联盟的格式更好。为了获得改进的数据集可读性，我们实际上需要合并三个表：FanGraphs、Baseball Reference和关键字查找。
from pybaseball import (cache, batting_stats_bref, batting_stats, playerid_reverse_lookup)import pandas as pdcache.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"\n相关性: {df[['OBP','SLG']].corr()['SLG']['OBP']:.3f}")
OBP SLGcount 362.000 362.000mean 0.320 0.415std 0.034 0.068min 0.234 0.22725% 0.300 0.36750% 0.318 0.41475% 0.340 0.460max 0.416 0.654相关性: 0.630让我们可视化这个联合分布，使用：
- 以国家联盟（NL）和美国联盟（AL）为颜色的球员散点图
- 球员的双变量核密度估计（KDE）图，使用高斯核对散点图进行平滑以估计密度
- 每个指标的边际KDE图
import matplotlib.pyplot as pltimport seaborn as snsg = 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("上垒率 vs. 长打率\n2023赛季，最少" f"{MIN_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()メソッドがありますが、例外やアルゴリズム間の引数の違いもあります。それぞれのアルゴリズムを個別に考慮しますが、各アルゴリズムをプレーヤーあたり（n = 2）の属性行列（m = 453）に適用します。そして、各プレーヤーだけでなく、各属性の範囲にまたがる値のグリッドも評価し、予測関数を視覚化するのに役立ちます。
決定境界を視覚化するためには、次の手順を実行する必要があります：
- 入力特徴値の2Dメッシュグリッドを作成する。
- 各点に
decision_functionを適用するために、グリッドのアンスタックが必要です。 - 予測をグリッドに再形成する。
- 予測をプロットする。
既存の観測範囲に加えて、パディングを含む200×200のグリッドを使用しますが、望む速度と解像度に合わせてグリッドを調整できます。
import numpy as npX = df[["OBP","SLG"]].to_numpy()GRID_RESOLUTION = 200disp_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.shapegrid_unstacked = np.c_[xx.ravel(), yy.ravel()]楕円包絡
楕円包絡の形状は、データの共分散行列によって決まります。共分散行列は、特徴iの分散を主対角線[i,i]に、特徴iと特徴jの共分散を[i,j]の位置に示します。共分散行列は外れ値に敏感なため、このアルゴリズムでは最小共分散行列決定法（MCD）推定量を使用します。これは、単峰性かつ対称分布に対して推奨されるもので、ランダム性は再現性のためのrandom_stateの入力によって決まります。この頑健な共分散行列は、後で再度役立つでしょう。
外れ値/内れ値のバイナリ分類ではなく、外れ値のスコアのランキングを比較したいため、decision_functionを使用してプレーヤーのスコアを計算します。
from sklearn.covariance import EllipticEnvelopeell = 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）に基づいています。各観測値からその最近傍までの総距離を計算し、局所密度を定義し、その後、各観測値の局所密度とその近傍の密度を比較します。近傍の密度よりもはるかに低い局所密度を持つ観測値は、外れ値と見なされます。
含める最近傍の数を選択する：KNNでは、規則として、K = sqrt(N)とし、Nは観測値の数です。この規則から、Kは20に近くなります（LOFのデフォルトKと偶然一致します）。過学習を減らすためにKを増やすか、適合不足を減らすためにKを減らすことができます。
K = int(np.sqrt(X.shape[0]))print(f"最近傍数K={K}を使用。")
最近傍数K=19を使用。距離尺度の選択：特徴量が相関し、異なる分散を持つため、ユークリッド距離はあまり意味がありません。特徴の尺度と相関を考慮したマハラノビス距離を使用します。
在计算马氏距离时，我们将使用鲁棒协方差矩阵。如果我们还没有通过椭圆包计算它，我们可以直接计算它。
from scipy.spatial.distance import pdist, squareform # 如果我们还没有椭圆包，我们可以计算鲁棒协方差：# from sklearn.covariance import MinCovDet# robust_cov = MinCovDet().fit(X).covariance_# 但我们可以直接从椭圆包中重复使用它：robust_cov = ell.covariance_print(f"鲁棒协方差矩阵：\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"马氏距离矩阵大小为 {D_mahal.shape}，例如：\n{np.round(D_mahal[:5,:5],3)}...\n...\n")
鲁棒协方差矩阵：[[0.00077 0.00095] [0.00095 0.00366]]马氏距离矩阵大小为(362, 362)，例如：[[0. 2.86 1.278 0.964 0.331] [2.86 0. 2.63 2.245 2.813] [1.278 2.63 0. 0.561 0.956] [0.964 2.245 0.561 0. 0.723] [0.331 2.813 0.956 0.723 0. ]]......
拟合局部异常因子：请注意，使用自定义距离矩阵要求我们向构造函数传递metric="precomputed"，然后将距离矩阵本身传递给fit方法。
还请注意，与其他算法不同，对于LOF，我们被禁止使用score_samples方法来对现有观测值进行评分；该方法只应用于新颖性检测。
from sklearn.neighbors import LocalOutlierFactorlof = LocalOutlierFactor(n_neighbors=K, metric="precomputed", novelty=True).fit(D_mahal)df["outlier_score_lof"] = lof.negative_outlier_factor_创建决策边界：因为我们使用了自定义距离度量，我们还必须计算网格中每个点到原始观测值的自定义距离。在之前，我们使用了空间度量pdist来计算单个集合中每个成员之间的成对距离，但现在我们使用cdist来返回第一个输入集合的每个成员到第二个输入集合的距离。
from scipy.spatial.distance import cdistD_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。我们将把经过缩放的输入数据输入到Nyström核逼近中，正如SGDOneClassSVM的文档所建议的。
from sklearn.pipeline import make_pipelinefrom sklearn.preprocessing import RobustScalerfrom sklearn.kernel_approximation import Nystroemfrom sklearn.linear_model import SGDOneClassSVMsuv = 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)孤立森林（Isolation Forest）
这种基于树的测算孤立性的方法执行随机递归划分。如果隔离给定观察所需的平均划分次数较低，则认为该观察是强有力的异常候选者。和随机森林以及其他基于树的模型一样，孤立森林不假设特征服从正态分布，也不要求特征进行缩放。默认情况下，它构建100棵树。我们的示例只使用了两个特征，因此我们没有启用特征抽样。
from sklearn.ensemble import IsolationForestiso = 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_textfrom sklearn.preprocessing import QuantileTransformerN_QUANTILES = 8 # 每个图表的颜色断点数N_CALLOUTS=15 # 每个图表中顶部异常值的标签数fig, axs = plt.subplots(2, 2, figsize=(12, 12), sharex=True, sharey=True)fig.suptitle("异常值识别算法比较",size=20)fig.supxlabel("出垒率（OBP）")fig.supylabel("垒打率（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], ["椭圆包络法","孤立森林", "局部离群因子","单类支持向量机"], 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) # 平均排名是任意的；我们只需要一个从大到小的倒序排列df["Rank_avg"] = df[model_rank_vars].mean(axis=1)print("倒序排列，最大的异常值...\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))正在倒计时最大的异常值…
倒计时到最大的异常值...
名称 AB PA H 2B 3B HR BB HBP SO OBP Pctl_OBP SLG Pctl_SLG 排名_ELL 排名_ISO 排名_LOF 排名_SVM Austin Barnes 178 200 32 5 0 2 17 2 43 0.256 2.6 0.242 0.6 19 7 25 12 J.D. Martinez 432 479 117 27 2 33 34 2 149 0.321 52.8 0.572 98.1 15 18 5 15 Yandy Diaz 525 600 173 35 0 22 65 8 94 0.410 99.2 0.522 95.4 13 15 13 10 Jose Siri 338 364 75 13 2 25 20 2 130 0.267 5.5 0.494 88.4 8 14 15 13 Juan Soto 568 708 156 32 1 35 132 2 129 0.410 99.2 0.519 95.0 12 13 11 11 Mookie Betts 584 693 179 40 1 39 96 8 107 0.408 98.6 0.579 98.3 7 10 20 7 Rob Refsnyder 202 243 50 9 1 1 33 5 47 0.365 90.5 0.317 6.6 5 19 2 14 Yordan Alvarez 410 496 120 24 1 31 69 13 92 0.407 98.3 0.583 98.6 6 9 18 6 Freddie Freeman 637 730 211 59 2 29 72 16 121 0.410 99.2 0.567 97.8 9 11 9 8 Matt Olson 608 720 172 27 3 54 104 4 167 0.389 96.5 0.604 99.2 11 6 7 9 Austin Hedges 185 212 34 5 0 1 11 2 47 0.234 0.3 0.227 0.3 10 1 4 3 Aaron Judge 367 458 98 16 0 37 88 0 130 0.406 98.1 0.613 99.4 3 5 6 4Ronald Acuna Jr. 643 735 217 35 4 41 80 9 84 0.416 100.0 0.596 98.9 2 3 10 2 Corey Seager 477 536 156 42 0 33 49 4 88 0.390 97.0 0.623 99.7 4 4 3 5 Shohei Ohtani 497 599 151 26 8 44 91 3 143 0.412 99.7 0.654 100.0 1 2 1 1分析和结论
看起来这四种实现在如何定义异常值方面基本是一致的，但在评分和易用性上存在一些明显的差异。
椭圆包络对于椭圆的短轴周围有更窄的轮廓，因此更倾向于突出那些与特征之间整体相关性相反的有趣选手。例如，光芒队的外野手何塞·西里（José Siri）在这个算法下更偏离正常，因为他的长打率（88百分位数）很高，而他的上垒率（5百分位数）很低，这与一个在边界球上用力挥棒，要么砸出好球，要么无法接触到球的果敢击球者一致。
椭圆包络也易于使用，无需配置，并提供稳健的协方差矩阵。如果你的数据是低维的，并且对其呈正态分布有合理的期望（通常并非如此），你可能想先尝试这种简单的方法。
单类支持向量机的轮廓更均匀间隔，因此更强调沿着整体相关方向的观察值。全明星一垒手弗雷迪·弗里曼（道奇队）和雅迪·迪亚兹（光芒队）在这个算法下的排名比其他算法更靠前，因为他们的长打率和上垒率都很出色（弗里曼分别是99百分位和97百分位，迪亚兹分别是99百分位和95百分位）。
使用RBF核函数时，需要额外进行标准化的步骤，但它在这个简单的例子中似乎表现良好，无需进一步调整。
局部异常因子注意到前面提到的“杰出球员群集”，它的轮廓呈现一个双峰的小轮廓（在图表上几乎看不到）。由于道奇队的外野手/二垒手穆基·贝茨（Mookie Betts）被其他出色的击球手包围，包括弗里曼、约旦·阿尔瓦雷斯和罗纳德·阿库尼亚·朱尼尔，他在LOF算法下的异常值排名仅为第20位，而在其他算法下的排名要靠前，至少在第10位。相反，勇士队的外野手马塞尔·奥祖纳（Marcell Ozuna）的长打率稍低，上垒率比贝茨低得多，但在LOF算法下更偏离正常，因为他的周围区域较为稀疏。
LOF是最耗时的实现，因为我们创建了用于拟合和评分的稳健距离矩阵。我们也可以花些时间调整K值。
隔离森林倾向于突出特征空间角落的观察值，因为切分点分布在各个特征上。替补捕手奥斯汀·海奇斯（Austin Hedges）在2023年为海盗队和游骑兵队效力，2024年加入守护者队，他在防守方面表现出色，但在长打率和上垒率两项指标中都是最差（至少200个打席）。海奇斯可以在上垒率或OPS中的一项指标上进行单一分割，这使他成为最强的异常值。隔离森林是唯一一个将大谷翔平视作最强异常值的算法：大谷翔平在上垒率上被罗纳德·阿库尼亚·朱尼尔略微压过，因此大谷翔平和阿库尼亚只能在一个特征上进行单一分割。
与常见的有监督树形学习器一样，隔离森林不进行外推，它更适合用于对受污染数据集进行异常值检测，而不是对无异常的数据集进行新颖性检测（在这种情况下，它不会将新的异常值评分比现有观察值更高）。
尽管隔离森林可以直接使用，但它未能将大谷翔平排为棒球界（甚至所有职业体育）最大的异常值，这说明任何异常值检测器的主要限制在于你用来拟合它的数据。
我们不仅省略了防守统计数据（对不起，奥斯汀·海奇斯），还没有包括投手的统计数据。因为投手现在都不再试图击球了…除了大谷翔平，他的赛季包括了最低的对手打击率（BAA）和棒球界第11的被自振率（ERA）（至少100局），一场完全比赛，以及一场他击出十个三振并且打了两个全垒得分的比赛。
有人认为大谷翔平是一个高级外星人冒充人类，但更有可能存在两个高级外星人冒充同一个人类。不幸的是，其中一个刚刚接受了肘部手术，将在2024年不会投球…但另一个刚刚签署了一份创纪录的10年7亿美元合同。多亏了异常值检测，现在我们可以知道为什么了！
