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

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}")

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()

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()]

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 = 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")

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)

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)

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)

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)
)