import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
import warnings
warnings.filterwarnings('ignore')

from scipy import stats
from sklearn.preprocessing import StandardScaler
from imblearn.over_sampling import SMOTE 
from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV, KFold
from sklearn.feature_selection import RFE 
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier, VotingClassifier, ExtraTreesClassifier
import xgboost as xgb
from sklearn.metrics import accuracy_score, auc, confusion_matrix, f1_score, \
    precision_score, recall_score, roc_curve  # 导入指标库
import prettytable
import sweetviz as sv # 自动eda
import toad 
from sklearn.model_selection import StratifiedKFold, cross_val_score  # 导入交叉检验算法

# 绘图初始化
%matplotlib inline
pd.set_option('display.max_columns', None) # 显示所有列
sns.set(style="ticks")
plt.rcParams['axes.unicode_minus']=False # 用来正常显示负号
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号

# 导入自定义模块
import sys
sys.path.append("/Users/heinrich/Desktop/Heinrich-blog/数据分析使用手册")
from keyIndicatorMapping import *

# 读取数据
raw_data = pd.read_csv('classification.csv')
raw_data.head()

# 缺失值填充，SMOTE方法限制非空
raw_data=raw_data.fillna(raw_data.mean())

# 数据集分割
X = raw_data[raw_data.columns.drop('churn')]
y = raw_data['churn']

# 标准化
scaler = StandardScaler() 
scale_data = scaler.fit_transform(X)  
X = pd.DataFrame(scale_data, columns = X.columns)

# 划分训练测试集
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=.3, random_state=0)

# 过采样
model_smote = SMOTE(random_state=0)  # 建立SMOTE模型对象
X_train, y_train = model_smote.fit_resample(X_train, y_train) 

%%time
# 初选分类模型
model_names = ['LR', 'SVC', 'RFC', 'XGBC']  # 不同模型的名称列表
model_lr = LogisticRegression(random_state=10) # 建立逻辑回归对象
model_svc = SVC(random_state=0, probability=True) # 建立支持向量机分类对象
model_rfc = RandomForestClassifier(random_state=10) # 建立随机森林分类对象
model_xgbc = xgb.XGBClassifier(use_label_encoder=False, eval_metric='auc', random_state=10) # 建立XGBC对象

# 模型拟合结果
model_list = [model_lr, model_svc, model_rfc, model_xgbc]  # 不同分类模型对象的集合
pre_y_list = [model.fit(X_train, y_train).predict(X_test) for model in model_list]  # 各个回归模型预测的y值列表

CPU times: user 2.49 s, sys: 125 ms, total: 2.62 s
Wall time: 843 ms

# 核心评估指标
metrics_dic = {
                'model_names':[],
                'auc':[],
                'ks':[],
                'accuracy':[],
                'precision':[],
                'recall':[],
                'f1':[]
            }
for model_name, model, pre_y in zip(model_names, model_list, pre_y_list):
    y_prob = model.predict_proba(X_test)  # 获得决策树的预测概率，返回各标签（即0，1）的概率
    fpr, tpr, thres = roc_curve(y_test, y_prob[:, 1])  # ROC y_score[:, 1]取标签为1的概率，这样画出来的roc曲线为正
    metrics_dic['model_names'].append(model_name)
    metrics_dic['auc'].append(auc(fpr, tpr))  # AUC
    metrics_dic['ks'].append(max(tpr - fpr)) # KS值
    metrics_dic['accuracy'].append(accuracy_score(y_test, pre_y))
    metrics_dic['precision'].append(precision_score(y_test, pre_y))
    metrics_dic['recall'].append(recall_score(y_test, pre_y))
    metrics_dic['f1'].append(f1_score(y_test, pre_y))
pd.DataFrame(metrics_dic)

%%time
# 建立组合评估器列表 均衡稳定性和准确性 这里只是演示，就将所有模型都纳入了
estimators = [('SVC', model_svc), ('RFC', model_rfc), ('XGBC', model_xgbc), ('LR', model_lr)]  
model_vot = VotingClassifier(estimators=estimators, voting='soft', weights=[1.1, 1.1, 0.9, 1.2],
                             n_jobs=-1)  # 建立组合评估模型
cv = StratifiedKFold(5)  # 设置交叉检验方法 分类算法常用交叉检验方法
cv_score = cross_val_score(model_vot, X_train, y_train, cv=cv, scoring='accuracy')  # 交叉检验
print('{:*^60}'.format('Cross val scores:'),'\n',cv_score) # 打印每次交叉检验得分
print('Mean scores is: %.2f' % cv_score.mean())  # 打印平均交叉检验得分

*********************Cross val scores:********************** 
 [0.73529412 0.7745098  0.85294118 0.85294118 0.87745098]
Mean scores is: 0.82
CPU times: user 2.38 s, sys: 432 ms, total: 2.81 s
Wall time: 5 s

# 模型训练
model_vot.fit(X_train, y_train)  # 模型训练

VotingClassifier(estimators=[('SVC', SVC(probability=True, random_state=0)),
                             ('RFC', RandomForestClassifier(random_state=10)),
                             ('XGBC',
                              XGBClassifier(base_score=0.5, booster='gbtree',
                                            colsample_bylevel=1,
                                            colsample_bynode=1,
                                            colsample_bytree=1,
                                            eval_metric='rmse', gamma=0,
                                            gpu_id=-1, importance_type='gain',
                                            interaction_constraints='',
                                            learning_rate=0.300000012,
                                            max...
                                            min_child_weight=1, missing=nan,
                                            monotone_constraints='()',
                                            n_estimators=100, n_jobs=8,
                                            num_parallel_tree=1,
                                            random_state=10, reg_alpha=0,
                                            reg_lambda=1, scale_pos_weight=1,
                                            subsample=1, tree_method='exact',
                                            use_label_encoder=False,
                                            validate_parameters=1,
                                            verbosity=None)),
                             ('LR', LogisticRegression(random_state=10))],
                 n_jobs=-1, voting='soft', weights=[1.1, 1.1, 0.9, 1.2])

model_confusion_metrics(model_vot, X_test, y_test, 'test')
model_core_metrics(model_vot, X_test, y_test, 'test')

confusion matrix for test
 +----------+--------------+--------------+
|          | prediction-0 | prediction-1 |
+----------+--------------+--------------+
| actual-0 |      53      |      31      |
| actual-1 |      37      |     179      |
+----------+--------------+--------------+
core metrics for test
 +-------+----------+-----------+--------+-------+-------+
|  auc  | accuracy | precision | recall |   f1  |   ks  |
+-------+----------+-----------+--------+-------+-------+
| 0.805 |  0.773   |   0.589   | 0.631  | 0.609 | 0.504 |
+-------+----------+-----------+--------+-------+-------+

# pip install shap

import shap   
# 初始化
shap.initjs()  

# 通过采样提高计算效率，但会导致准确率降低。表现在base_value与mean(model.predict_proba(X))存在差异，不建议K太小
# X_test_summary = shap.sample(X_test, 200)
# X_test_summary = shap.kmeans(X_test, 150)
explainer = shap.KernelExplainer(model_vot.predict_proba, X_test)
shap_values = explainer.shap_values(X_test, nsamples = 10)

Using 300 background data samples could cause slower run times. Consider using shap.sample(data, K) or shap.kmeans(data, K) to summarize the background as K samples.

# 单样本查看-1概率较高的样本 # 208
shap.force_plot(base_value=explainer.expected_value[1],
                shap_values=shap_values[1][208],
                features = X_test.iloc[208,:]
               )

# 验证base_value
print('所有样本预测标签1的概率均值：',model_vot.predict_proba(X_test)[:,1].mean())
print('base_value:',explainer.expected_value[1])

所有样本预测标签1的概率均值：0.3519852365700774
base_value: 0.35198523657007774

# 验证单一样本
i=208
fx=model_vot.predict_proba(X_test)[:,1][i]
da=X_test.iloc[i,:]
sv=fx-explainer.expected_value[1]
sv_val=shap_values[1][i].sum()
print('f(x):',fx)
print('shap_values:',sv,sv_val)

f(x): 0.9264517406651224
shap_values: 0.5744665040950446 0.5744665040950446

# 特征重要程度
shap.summary_plot(shap_values[1],X_test,max_display=10,plot_type="bar")

# 特征与样本蜂窝图
shap.summary_plot(shap_values[1],X_test,max_display=10)

# 单特征预测结果
shap.dependence_plot("retention_days", shap_values[1], X_test, interaction_index=None)

# 双特征交叉影响
shap.dependence_plot("retention_days", shap_values[1], X_test, interaction_index='level')