特征漂移指标 PSI

size = 5000
# 期望数据
p2 = np.random.normal(loc = 3, scale = 1, size = size)
# 实际数据
a2 = np.random.normal(loc = 3.5, scale = 0.75, size = size)

num_bins = 10
eps = 1e-4

min_val = min(min(p2), min(a2))
max_val = max(max(a2), max(p2))
bins = [min_val + (max_val - min_val)*(i)/num_bins for i in range(num_bins+1)]
bins[0] = min_val - eps # 修正下界
bins[-1] = max_val + eps # 修正上界

print(bins)
# [-0.4810252475657688, 0.229173950184835, 0.9392731479354388, 1.649372345686043, 2.3594715434366464, 3.0695707411872504, 3.7796699389378547, 4.489769136688458, 5.199868334439062, 5.909967532189666, 6.620166729940269]

# 计算数组总的元素属于哪个分箱
bins_p2 = pd.cut(p2, bins = bins, labels = range(1,num_bins+1))
# 将元素和封箱号对齐
df_p2 = pd.DataFrame({'p2': p2, 'bin': bins_p2})
# 统计每个封箱中的元素数量
grp_p2 = df_p2.groupby('bin').count()
# 计算每个分箱中元素数量占总数的百分比
grp_p2['percent_p2'] = grp_p2['p2'] / sum(grp_p2['p2'])

# 对实际数据也进行分享操作
bins_a2 = pd.cut(a2, bins = bins, labels = range(1,num_bins+1))
df_a2 = pd.DataFrame({'a2': a2, 'bin': bins_a2})
grp_a2 = df_a2.groupby('bin').count()
grp_a2['percent_a2'] = grp_a2['a2'] / sum(grp_a2['a2'])

# 比较 p2 和 a2 的分箱数据
psi_df = grp_p2.join(grp_a2, on = "bin", how = "inner")
print(psi_df)

# 当 percent_p* 是 0 时，给其加上 eps。防止 0 参与运算抛异常。
psi_df['percent_p2'] = psi_df['percent_p2'].apply(lambda x: eps if x == 0 else x)
psi_df['percent_a2'] = psi_df['percent_a2'].apply(lambda x: eps if x == 0 else x)
# 计算每行的 psi
psi_df['psi'] = (psi_df['percent_p2'] - psi_df['percent_a2']) * np.log(psi_df['percent_p2'] / psi_df['percent_a2'])

print(psi_df)

psi = psi_df['psi'].sum()

print(psi)
# 0.4533650280982507

size = 5000
# 合成训练需要的数据
x1 = np.random.normal(loc = 0, scale = 2, size = size)
x2 = np.random.normal(loc = 3, scale = 1, size = size)
x3 = np.random.normal(loc = 5, scale = 2, size = size)
y = [-2*x1 + 3.1415*x2 + 2.7183*x3 + 1.6180 for (x1, x2, x3) in zip(x1,x2,x3)]

# 划分训练集和验证集
X = pd.DataFrame({'x1': x1, 'x2': x2, 'x3': x3})
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2)

# 训练模型（该模型尽可能拟合出 y = a1*x1 + a2*x2 + a3*x3 + b 表达式）
model = ElasticNet()
model.fit(X_train, y_train)

# 验证结果（模型线下的评测结果）
y_pred = model.predict(X_test)

# 再合成生产环境（线上）数据
z1 = np.random.normal(loc = 0, scale = 2, size = size)  
z2 = np.random.normal(loc = 3.5, scale = 0.75, size = size)  
z3 = np.random.normal(loc = 8, scale = 3, size = size)
# 计算模型的线上推理结果
Z = pd.DataFrame({'x1': z1, 'x2': z2, 'x3': z3})
z_pred = model.predict(Z)

"""
此时线下模型效果 y_pred 和线上模型效果 z_pred 都是一维的，令 p2 = y_pred, a2 = z_pred。
套用“PSI 代码实现部分”就可以计算出群体性稳定指标 PSI，最终判定出模型稳定性处于什么程度。
"""

令 p2 = x1, a2 = z1, 计算出 psi1;
令 p2 = x2, a2 = z2, 计算出 psi2;
令 p2 = x3, a2 = z3, 计算出 psi3;
通过 psi* 就能看出是哪个特征引起的不稳定。