基於最優分箱的psi值計算-python實現

一、需求：

篩選目標相關變數時，考慮變數分佈在時間變化上的穩定性；

二、解決思路：

將生產數據按照時間切割
計算訓練集合var變數最優分箱
計算各個數據塊在最優分箱的分佈頻率
計算不同數據塊上變數var的psi值。

三、使用數據：

train集合：包含檢驗變數var、目標變數target；

p集合（生產數據）：包含檢驗變數var、樣本觀測時間date；

四、設計函數：

1、函數1：使用訓練集合對觀測變數var進行最優分箱

輸入：train集合（dataframe格式）

輸出：最優分箱閾值（list格式）

2、函數2：對p集合按照date倒序切割，並存儲到pp列表裡

輸入：p集合（dataframe格式）

輸出：pp（list格式，）

3、函數3：按照閾值統計分佈頻率

輸入：閾值

輸出：閾值分佈頻率

4、函數4：對兩個長度相等的list列表，計算psi值

輸入：原分佈頻率base、比較分佈頻率p_list

輸出：psi值

5、主函數：調用函數1~4，完成psi值的統計

輸入：訓練集、模型目標、生產數據、分組欄位、觀察

輸出：橫軸為分組縱軸為觀察欄位的表格

五、代碼部分：

#卡方分箱 def chi3(arr): assert(arr.ndim==2) R_N = arr.sum(axis=1) C_N = arr.sum(axis=0) N = arr.sum() E = np.ones(arr.shape)* C_N / N E = (E.T * R_N).T square = (arr-E)**2 / E square[E==0] = 0 v = square.sum() return v def chiMerge(df,col,target,max_groups=None,threshold=None): freq_tab = pd.crosstab(df[col],df[target]) freq = freq_tab.values cutoffs = freq_tab.index.values if max_groups is None: if threshold is None: cls_num = freq.shape[-1] threshold = chi2.isf(0.05,df= cls_num - 1)

while True:
minvalue = None
minidx = None
for i in range(len(freq) - 1):
v = chi3(freq[i:i+2])
if minvalue is None or (minvalue > v): #小於當前最小卡方，更新最小值
minvalue = v
minidx = i
if (max_groups is not None and max_groups< len(freq) ) or (threshold is not None and minvalue < threshold):
tmp = freq[minidx] + freq[minidx+1]
freq[minidx] = tmp
freq = np.delete(freq,minidx+1,0)
cutoffs = np.delete(cutoffs,minidx+1,0)
else: #最小卡方值不小於閾值，停止合併。
break
return cutoffs

#切割陪跑數據
def get_plist(p,by):
p.sort_values(by=[by],ascending=False,inplace=True)
rk=list(p[by].unique())
pp=[]
for i in range(len(rk)):
pp.append(p[p[by]==rk[i]])
return pp,rk

#按照閾值統計分佈頻率：
def get_percent(d,n,cutoffs): #計算d集合n特徵在最優分箱上的分佈
s=[]
p=[]
for m in range(len(cutoffs)):
if m==(len(cutoffs)-1):
s.append(d[d[n]>=cutoffs[m]][n].count())
else:
s.append(d[d[n]>=cutoffs[m]][d[n]<cutoffs[m+1]][n].count())
for i in s:
p.append(i/sum(s))
return p

def psi(s1,s2): #以s1概率分佈為標準數據，計算s2概率分佈的psi值
psi=0
for i in range(len(s1)):
if s2[i]==0:
s2[i]=0.000001
if s1[i]==0:
s1[i]=0.000001
p=((s2[i]-s1[i])*(math.log((s2[i]/s1[i]))))
psi=psi+p
return psi

def main(train,target,p,by,columns):
pp,by=get_plist(p,by)
result=pd.DataFrame(index=by,columns=columns)
for n in columns:
print(n)
cutoffs=chiMerge(train,n,target,max_groups=8)
for d in range(len(pp)):
base_p=get_percent(pp[1],n,cutoffs) #使用第二部分作為比較對象（時間最近的完整數據）
d_p=get_percent(pp[d],n,cutoffs)
psi1=psi(base_p,d_p)
print(by[d],psi1)
result.loc[by[d],n]=psi1
return result
#示例：
psi=main(X1,TARGET,p,DATE,n)
print(psi)