譯者 | 朱先忠
審校 | 重樓
異常值檢測是一種無監督的機器學習任務，用于識別給定數據集中的異常（即“異常觀測”）。在大量現實世界中，當我們的可用數據集已經被異?！拔廴尽睍r，異常值檢測任務對于整個機器學習環節來說是非常有幫助的。當前，開源框架Scikit learn已經實現了幾種異常值檢測算法，在我們持有未受污染的基準數據的情況下，我們也可以使用這些算法進行新穎檢測。這是一項半監督任務，可以預測新的觀測結果是否為異常值。
概述
本文中，我們將比較的四種異常值檢測算法如下：
- 橢圓包絡（Elliptic Envelope）適用于低維正態分布數據。顧名思義，它使用多元正態分布來創建一個距離度量指標，以便將異常值與內部值區分開來。
- 局部異常因子（Local Outlier Factor）是一種針對觀測的局部密度與其鄰居的局部密度的比較。密度比其鄰居低得多的觀測被視為異常值。
- 具有隨機梯度下降的單類支持向量機（One-Class Support Vector Machine with Stochastic Gradient Descent）是一類SVM的O（n）近似解。請注意，O（n2）數量級的單類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）中。然后，我們不僅會為每個球員打分，還會為每個屬性范圍內的值網格打分，以幫助我們可視化預測函數。
為了可視化決策邊界，我們需要采取以下步驟：
- 創建輸入特征值的二維網格。
- 將decision_function應用于網格上的每個點，這需要拆開網格。
- 將預測重新成形為一個網格。
- 繪制預測。
在此，我們將使用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）算法。我們計算每個觀測到其最近鄰居的總距離，以定義局部密度，然后將每個觀測的局部密度與其鄰近的局部密度進行比較。局部密度遠小于其鄰近的觀測結果即被視為異常值。
選擇要包括的鄰近數量：在KNN中，經驗法則是讓K=sqrt（N），其中N是您的觀察計數。根據這個規則，我們得到了一個接近20的K（這恰好是LOF的默認K）。您可以分別增加或減少K以減少過擬合或欠擬合。
K = int(np.sqrt(X.shape[0]))
print(f"Using K={K} nearest neighbors.")
選擇一個距離指標：請注意，我們的特征是相關的，并且具有不同的方差，因此歐幾里得距離不是很有意義。我們將使用Mahalanobis距離，這種距離更有助于體現特征比例和相關性。
在計算馬氏距離時，我們將使用魯棒性強的協方差矩陣。如果我們還沒有通過橢圓包絡計算它，那么我們可以直接計算它。
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")
擬合局部異常因子：請注意，使用自定義距離矩陣需要將metric=“precomputed”傳遞給構造函數，然后將距離矩陣本身傳遞給fit方法。（有關更多詳細信息，請參閱官方文檔：https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.LocalOutlierFactor.html#sklearn.neighbors.LocalOutlierFactor.fit）
還要注意，與其他算法不同，使用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_創建決策邊界：因為我們使用了自定義距離度量，所以我們還必須計算網格中每個點與原始觀測值之間的自定義距離。以前，我們使用空間測度pdist來表示單個集合的每個成員之間的成對距離，但現在我們使用cdist來返回從第一組輸入的每個成員到第二組輸入的各個成員的距離。
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)
)
分析與結論
看起來上面四種實現方案在如何定義異常值方面基本保持一致，但在分值和易用性方面存在一些明顯差異。
橢圓包絡方法：在橢圓的短軸周圍有更窄的輪廓，因此它傾向于突出那些與特征之間的整體相關性相反的有趣球員。例如，在這種算法下，光芒隊的外野手JoséSiri更像是一個異類，因為他的SLG很高（第88百分位），而OBP很低（僅為第5百分位）。
橢圓包絡在沒有配置的情況下也很容易使用，并且它提供了具有魯棒性的協方差矩陣。如果您有低維數據，并且合理地期望它是正態分布的（通常情況并非如此），那么您可能想先嘗試這種簡單的方法。
單分類算法：具有更均勻分布的輪廓，因此它傾向于比橢圓包絡更強調沿整體相關方向的觀測。全明星一壘手Freddie Freeman（道奇隊）和Yandy Diaz（光芒隊）在該算法下的排名比其他算法下的更靠前，因為他們的SLG和OBP都很出色（Freeman分別為第99和97百分位，Diaz分別為第98和95百分位）。
值得注意的是，RBF內核需要一個額外的步驟來進行標準化，但在這個簡單的例子中，它似乎也能很好地工作，而無需進行微調。
局部異常因子：在前面提到的“優秀集群”上出現了一個小的雙峰輪廓（在圖表中幾乎看不到）。由于道奇隊的外野手/二壘手穆基·貝茨周圍都是其他優秀的擊球手，包括弗里曼、約丹·阿爾瓦雷斯和小羅納德·阿庫尼亞，他在LOF下僅排名第20，而在其他算法下排名第10或更強。相反，勇士隊外野手Marcel Ozuna的SLG和OBP略低于Betts；但在LOF下，他更像是一個異類，因為他的鄰近密度較小。
另外，值得注意的是，LOF實現是最耗時的，因為我們要創建的是用于擬合和評分的穩健的距離矩陣。當然，我們本可以花一些時間來調整K。
孤立森林：傾向于強調在特征空間的角落進行觀察，因為分割分布在特征之間。替補捕手奧斯汀·赫奇斯于2023年效力于海盜隊和流浪者隊，并于2024年與守護者隊簽約。他具有很強的防守能力，但在SLG和OBP中都是最差的擊球手（至少有200次板上擊球）。赫奇斯可以在OBP或OPS上單獨分離，使他成為最強的異類。隔離森林是唯一一種沒有將大谷昭平列為最強異常值的算法：由于大谷昭平被小羅納德·阿庫尼亞在OBP中淘汰；因此，大谷昭和阿庫尼亞只能在一個特征上進行單獨分離。
與常見的基于監督樹的學習器一樣，隔離森林不進行外推分析，這使得它更適合于擬合受污染的數據集進行異常值檢測，而不是擬合無異常的數據集用于新穎性檢測（在這種情況下，它不會比現有觀察更有力地對新的異常值進行評分）。
盡管“隔離森林”開箱即用效果很好，但它未能將大谷昭平列為棒球（可能是所有職業運動）中最大的異常值，這說明了任何異常值檢測器的主要局限性：用于擬合它的數據。
需要說明的是，在我們的實驗數據中，我們不僅省略了防守數據（對不起，奧斯汀·赫奇斯），也沒有考慮投球數據。因為投手們甚至不再嘗試擊球了……除了Ohtani，他的賽季包括棒球史上第二好的打擊率（BAA）和第11好的平均得分（ERA）（至少投了100局），一場完整的比賽被淘汰，以及一場他三振十名擊球手并擊出兩支全壘打的比賽。
有人認為大谷昭平是一個模仿人類的高級外星人，但似乎更有可能有兩個模仿同一個人的高級外星人。不幸的是，其中一人剛剛做了肘部手術，2024年不會投球了……但另一人剛剛簽署了一份創紀錄的10年7億美元的合同。哈哈，得益于異常值檢測，現在我們終于可以明白其中有關的原因了！
譯者介紹
朱先忠，51CTO社區編輯，51CTO專家博客、講師，濰坊一所高校計算機教師，自由編程界老兵一枚。
原文標題：Comparing Outlier Detection Methods，作者：John Andrews
