「アウトライア検出手法の比較」
「アウトライア検出手法の比較に関する研究」
2023年のメジャーリーグベースボールの打撃統計を使用する
異常検出は、与えられたデータセット内の異常（異常な観測）を特定するための教師なしの機械学習タスクです。このタスクは、利用可能なデータセットが既に「汚染」されている現実のケースで役立ちます。Scikit-learn では、いくつかの異常検出アルゴリズムが実装されており、汚染されていないベースラインがある場合には、これらのアルゴリズムを新規性検出のためにも使用することができます。新規性検出は半教師ありのタスクであり、新たな観測が異常であるかどうかを予測します。
概要
比較する4つの異常検出アルゴリズムは次のとおりです:
- 楕円包絡線法（Elliptic Envelope）は、次元数が低い正規分布データに適しています。その名前が示すように、多変量正規分布を使用して外れ値と内れ値を分離するための距離尺度を作成します。
- 局所外れ値係数（Local Outlier Factor）は、観測値の局所密度とその近傍の密度を比較します。近傍よりもはるかに密度が低い観測値は外れ値と見なされます。
- 確率的勾配降下法付きワンクラスサポートベクターマシン（One-Class Support Vector Machine (SVM) with Stochastic Gradient Descent (SGD)）はワンクラスSVMのO(n)近似解です。なお、O(n²)のワンクラスSVMは小さなデータセットの例ではうまく動作しますが、実際の使用例では実用的ではありません。
- 孤立した木（Isolation Forest）は、ランダムな分割によって外れ値をより速く分離する木ベースの手法です。
私たちのタスクは教師なしのため、これらのアルゴリズムの正確さを比較するための正解データはありません。代わりに、それぞれの結果（特に選手の順位）がどのように異なるかを見て、彼らの振る舞いや制限についての直感を得ることで、いつどのアルゴリズムを優先すべきかを知ることが目的です。
2023年のメジャーリーグベースボール（MLB）シーズンのバッターのパフォーマンスから、いくつかの技術を比較してみましょう。
- 出塁率（OBP）：打者がベースに到達する割合（ヒット、四球、死球を含む）
- 長打率（SLG）：打席あたりの総塁数の平均
打者のパフォーマンスのさまざまな洗練された指標がありますが、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、およびキールックアップの3つのテーブルを結合する必要があります。
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":"リーグ", "Tm":"チーム"} )df_bref["リーグ"] = \ df_bref["リーグ"].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="チーム" ).merge(key_mapping, how="inner", on="IDfg" ).merge(df_bref[["mlbID","リーグ","チーム"]], how="inner", on="mlbID" ).sort_values(["リーグ","チーム","名前"])データの探索
まず、これらの指標は平均と分散が異なり、相関があります。また、各指標は比較的対称的であり、中央値は平均値に近いです。
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="リーグ", 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["名前"], (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
meshgridを作成します。 meshgrid上の各点に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は観測値の数です。このルールから、約20に近いKを得ることができます（ちょうどLOFのデフォルトのKとなります）。過学習を減らすためにKを増やしたり、適合不足を減らすためにKを減らすことができます。
K = int(np.sqrt(X.shape[0]))print(f"Using K={K} nearest neighbors.")
Using K=19 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"頑健な共分散行列:\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}のマハラノビス距離行列、" f"例えば:\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. ]]......LOFの適合: カスタム距離行列を使用するためには、コンストラクタに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_決定境界の作成: カスタム距離メトリックを使用したため、グリッド内の各点から元の観測値までのカスタム距離を計算する必要があります。以前は単一のセットの各メンバー間のペアワイズ距離のために空間計量を使用していましたが、今回は第一セットの各メンバーから第二セットの各メンバーまでの距離を返すために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)孤立森林
このツリーベースのアプローチは、ランダムな再帰分割を行い、孤立を測定します。特定の観測値を孤立させるために必要な分割の平均数が低い場合、その観測値はより強力な外れ候補と見なされます。ランダムフォレストや他のツリーベースのモデルと同様に、孤立森林は特徴が正規分布していることを仮定せず、特徴をスケーリングする必要もありません。デフォルトでは、100本の木を作成します。この例では2つの特徴のみを使用しているため、特徴のサンプリングは有効にしていません。
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("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) # ランクの平均化は任意です。ただし、カウントダウンの順序が必要です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 Rank_ELL Rank_ISO Rank_LOF Rank_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分析と結論
4つの実装は、外れ値の定義方法についてほぼ同意しているようですが、スコアや使いやすさには顕著な違いがあります。
楕円エンベロープは楕円の短軸周りに狭い輪郭を持っているため、特徴間の全体的な相関に反する興味深いプレーヤーを強調する傾向があります。例えば、レイズの外野手、ジョゼ・シリは、このアルゴリズムでは外れ値として扱われます。彼のSLG（88パーセンタイル）が低いOBP（5パーセンタイル）とは、境界線の球に積極的に振り、それを打ち払うかあまり力のない打撃をする攻撃的な打者と一致しています。
楕円エンベロープは設定なしで使いやすく、頑健な共分散行列を提供します。低次元のデータで正規分布されることを合理的に期待できる場合（通常はそうではない），このシンプルな手法を最初に試してみる価値があります。
ワンクラスSVMはより均一に配置された輪郭を持ち、楕円エンベロープよりも全体的な相関方向の観察に重点を置く傾向があります。オールスターの一塁手、フレディ・フリーマン（ドジャース）とヤンディ・ディアズ（レイズ）は、このアルゴリズムのもとで他のアルゴリズムよりも高くランクされます。彼らのSLGとOBPはどちらも優れています（フリーマンはそれぞれ99パーセンタイルと97パーセンタイル，ディアズはそれぞれ99パーセンタイルと95パーセンタイル）。
RBFカーネルは標準化のために追加の手順が必要でしたが、微調整なしでもこのシンプルな例ではうまく動作したようです。
ローカル外れ値ファクタは、前述の「卓越のクラスタ」を検出しました（グラフではほとんど目に見えません）。ドジャースの外野手／セカンドベースマンであるムーキー・ベッツは、フリーマン、ヨルダン・アルバレス、ロナルド・アクーニャ・ジュニアなど優れた打者たちに囲まれているため、他のアルゴリズムでは10番目以上の外れ値となりますが、LOFでは20番目に強い外れ値としてランクされます。逆に、ブレーブスの外野手であるマーセル・オズナは、ベッツよりもわずかに低いSLGとかなり低いOBPを持っていますが、彼の近隣はより密度が低いため、LOFではより外れ値となります。
LOFは、適合とスコアリングのために頑健な距離行列を作成したため、実装に最も時間がかかりました。Kの微調整にも時間を使うことができました。
孤立した木は特徴空間の角に観測を重点的に置く傾向があります。2023年にピレートスとレンジャーズでプレーし、2024年にガーディアンズと契約したバックアップキャッチャーのオースティン・ヘッジスは、守備では強いがSLGとOBPの両方で最低（少なくとも200打席）です。ヘッジスはOBPまたはOPSのいずれかの一つで単一の分割で孤立できるため、最も強い外れ値として扱われます。孤立した木は唯一、最も強い外れ値として大谷翔平をランク付けしなかったアルゴリズムです。大谷翔平はOBPでロナルド・アクーニャ・ジュニアに抜かれたため、両者とも単一の特徴量の分割で孤立できます。
一般的な教師付き木ベースの学習アルゴリズムと同様に、孤立した木は外挿を行わないため、新規データや既存の観測よりも新しい外れ値を強くスコアリングすることに適しています。
孤立した木は、そのまま使用してうまく機能しましたが、野球では（おそらくすべてのプロスポーツでも）最も強力な外れ値としての大谷翔平をランク付けできなかったことは、任意の外れ値検出器の主な制限を示しています：それを適合させるために使用するデータ。
防御統計を省略しただけでなく（すみません、オースティン・ヘッジス）、投手の統計も含めていませんでした。なぜなら、ピッチャーはもはや打つことを試みないからです…ただし、大谷翔平は例外です。彼のシーズンは他の選手よりも優れた打率対戦（BAA）と11位の防御率（ERA）（最低100イニング）を含み、完封勝利もあり、10人を三振に取り、2本のホームランも放ちました。
大谷翔平は人間をなりすましている高度な地球外生命体であるという指摘もありますが、同じ人間をなりすましている高度な地球外生命体が2人いる可能性の方が高いようです。残念ながら、そのうちの1人は肘の手術を受け、2024年に投げることはありません…しかし、もう1人は新記録となる10年間で7億ドルの契約を結びました。そして、外れ値検出によって、なぜ彼がそうなるのかがわかるようになりました！
We will continue to update VoAGI; if you have any questions or suggestions, please contact us!
Was this article helpful?
93 out of 132 found this helpful
Related articles
