图 4 - KL散度具有非对称性
因此，KL散度是单向描述信息熵差异。为了加强大家理解，接下来以数值案例介绍相对熵如何计算。
假如一个字符发射器，随机发出0和1两种字符，真实发出概率分布为A，但实际不知道A的具体分布。通过观察，得到概率分布B与C，各个分布的具体情况如下：
图 5 - 三种数据分布
可以计算出得到KL散度如下：
(3)KL(A||B)=1/2log(1/21/4) 1/2log(1/23/4=1/2log(4/3)KL(A||C)=1/2log(1/21/8) 1/2log(1/27/8=1/2log(16/7)
由上式可知，相对熵KL(A||C) > KL(A||B)，说明A和B之间的概率分布在信息量角度更为接近。而通过概率分布可视化观察，我们也认为A和B更为接近，两者吻合。
Part 4. 相对熵与PSI之间的关系
接下来，我们从数学上来分析相对熵和PSI之间的关系。
(4)psi=∑i=1n(Ai−Ei)∗ln(Ai/Ei)psi=∑i=1nAi∗ln(Ai/Ei) ∑i=1nEi∗ln(Ei/Ai)
将PSI计算公式变形后可以分解为2项，其中：
第1项：实际分布（A）与预期分布（E）之间的KL散度—— KL(A||E)
第2项：预期分布（E）与实际分布（A）之间的KL散度—— KL(E||A)
因此，PSI本质上是实际分布（A）与预期分布（E）的KL散度的一个对称化操作。其双向计算相对熵，并把两部分相对熵相加，从而更为全面地描述两个分布的差异。
至此，已经回答了上文的问题Q1——为什么要在计算公式中引入对数？
可能有人还会问，那为什么不把PSI定义成下式呢？
(5)KL(A||E) KL(E||A)2
按笔者理解，一是为了保证数学公式的优雅；二是只要在同一尺度上比较，同时放大2倍也无所谓。
再回到上文的问题Q2——为什么log的底数是自然常数e，而不是其他常数（比如2）？按笔者理解，两者都可以，从信息量角度而言，底数为2或许更为合理。
至于问题Q3——为什么不定义为“SUM(实际占比 - 预期占比）”？这是因为有正有负，会存在相互抵消的现象。另一方面，也确实没有从相对熵角度来定义更为合理。
图 6 - PSI与KL散度之间的关系
Part 5. PSI指标的业务应用
我们在业务上一般怎么用呢？一般以训练集（INS）的样本分布作为预期分布，进而跨时间窗按月/周来计算PSI，得到如图6所示的Monthly PSI Report，进而剔除不稳定的变量。同理，在模型上线部署后，也将通过PSI曲线报表来观察模型的稳定性。
入模变量保证稳定性，变量监控
模型分数保证稳定性，模型监控
根据建模经验，给出一些建议：
1. 实际评估需要分不同粒度：时间粒度（按月、按样本集）、订单层次（放贷层、申请层）、人群（若没有分群建模，可忽略）。
2. 先在放贷样本上计算PSI，剔除不稳定的特征；再对申请样本抽样（可能数据太大），计算PSI再次筛选。之前犯的错误就是只在放贷样本上评估，后来在全量申请订单上评估时发现并不稳定，导致返工。
3. 时间窗尽可能至今为止，有可能建模时间窗稳定，但近期时间窗出现不稳定。
4. PSI只是一个宏观的指标，建议先看变量数据分布（EDD），看分位数跨时间变化来检验数据质量。我们无法得知PSI上升时，数据分布是左偏还是右偏。因此，建议把PSI计算细节也予以保留，便于在模型不稳定时，第一时间排查问题。
图 7 - Monthly PSI Report
Part 6. PSI的计算代码（Python）
import math
import numpy as np
import pandas as pd
def calculate_psi(base_list, test_list, bins=20, min_sample=10):
try:
base_df = pd.DataFrame(base_list, columns=['score'])
test_df = pd.DataFrame(test_list, columns=['score'])
# 1.去除缺失值后，统计两个分布的样本量
base_notnull_cnt = len(list(base_df['score'].dropna()))
test_notnull_cnt = len(list(test_df['score'].dropna()))
# 空分箱
base_null_cnt = len(base_df) - base_notnull_cnt
test_null_cnt = len(test_df) - test_notnull_cnt
# 2.最小分箱数
q_list = []
if type(bins) == int:
bin_num = min(bins, int(base_notnull_cnt / min_sample))
q_list = [x / bin_num for x in range(1, bin_num)]
break_list = []
for q in q_list:
bk = base_df['score'].quantile(q)
break_list.append(bk)
break_list = sorted(list(set(break_list))) # 去重复后排序
score_bin_list = [-np.inf] break_list [np.inf]
else:
score_bin_list = bins
# 4.统计各分箱内的样本量
base_cnt_list = [base_null_cnt]
test_cnt_list = [test_null_cnt]
bucket_list = ["MISSING"]
for i in range(len(score_bin_list)-1):
left = round(score_bin_list[i 0], 4)
right = round(score_bin_list[i 1], 4)
bucket_list.append("(" str(left) ',' str(right) ']')
base_cnt = base_df[(base_df.score > left) & (base_df.score <= right)].shape[0]
base_cnt_list.append(base_cnt)
test_cnt = test_df[(test_df.score > left) & (test_df.score <= right)].shape[0]
test_cnt_list.append(test_cnt)
# 5.汇总统计结果
stat_df = pd.DataFrame({"bucket": bucket_list, "base_cnt": base_cnt_list, "test_cnt": test_cnt_list})
stat_df['base_dist'] = stat_df['base_cnt'] / len(base_df)
stat_df['test_dist'] = stat_df['test_cnt'] / len(test_df)
def sub_psi(row):
# 6.计算PSI
base_list = row['base_dist']
test_dist = row['test_dist']
# 处理某分箱内样本量为0的情况
if base_list == 0 and test_dist == 0:
return 0
elif base_list == 0 and test_dist > 0:
base_list = 1 / base_notnull_cnt
elif base_list > 0 and test_dist == 0:
test_dist = 1 / test_notnull_cnt
return (test_dist - base_list) * np.log(test_dist / base_list)
stat_df['psi'] = stat_df.apply(lambda row: sub_psi(row), axis=1)
stat_df = stat_df[['bucket', 'base_cnt', 'base_dist', 'test_cnt', 'test_dist', 'psi']]
psi = stat_df['psi'].sum()
except:
print('error!!!')
psi = np.nan
stat_df = None
return psi, stat_df
运行演示：
>>> psi, stat_df = calculate_psi(base_list=list(df[df['draw_month'] == '2020-05'][var]),
test_list=list(df[df['draw_month'] == '2021-02'][var]),
bins=20, min_sample=10)
>>> stat_df
bucket base_cnt base_dist test_cnt test_dist psi
0 MISSING 0 0.000000 0 0.000000 0.000000e 00
1 (-inf,0.0199] 238 0.046777 511 0.065454 6.274898e-03
2 (0.0199,0.0217] 266 0.052280 477 0.061099 1.374765e-03
3 (0.0217,0.0235] 249 0.048939 467 0.059818 2.183942e-03
4 (0.0235,0.0249] 262 0.051494 380 0.048674 1.587606e-04
5 (0.0249,0.0265] 280 0.055031 453 0.058025 1.585503e-04
6 (0.0265,0.0279] 230 0.045204 329 0.042142 2.148738e-04
7 (0.0279,0.0297] 251 0.049332 410 0.052517 1.992932e-04
8 (0.0297,0.0316] 254 0.049921 388 0.049699 9.929669e-07
9 (0.0316,0.0336] 248 0.048742 386 0.049443 1.000043e-05
10 (0.0336,0.0358] 257 0.050511 395 0.050596 1.416201e-07
11 (0.0358,0.0382] 263 0.051690 363 0.046497 5.499267e-04
12 (0.0382,0.0405] 249 0.048939 343 0.043935 5.396963e-04
13 (0.0405,0.0431] 258 0.050708 316 0.040476 2.305601e-03
14 (0.0431,0.046] 256 0.050314 319 0.040861 1.967526e-03
15 (0.046,0.0492] 255 0.050118 321 0.041117 1.781819e-03
16 (0.0492,0.0528] 252 0.049528 314 0.040220 1.937663e-03
17 (0.0528,0.0578] 255 0.050118 391 0.050083 2.398626e-08
18 (0.0578,0.0641] 262 0.051494 369 0.047265 3.623085e-04
19 (0.0641,0.0749] 249 0.048939 452 0.057897 1.505794e-03
20 (0.0749,inf] 254 0.049921 423 0.054182 3.489647e-04
加强对各类指标的理解，才能巩固风控模型基本功。
来源：是汪在路上（知乎）
