01 集成学习 - 概述、Bagging - 随机森林、袋外错误率
02 集成学习 - 特征重要度、Extra Tree、TRTE、IForest、随机森林总结
03 集成学习 - Boosting - AdaBoost算法原理
04 集成学习 - Boosting - AdaBoost算法构建
05 集成学习 - Boosting - GBDT初探
06 集成学习 - Boosting - GBDT算法原理、总结

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn import tree
# 引入了集成学习的随机森林库
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
from sklearn.preprocessing import Imputer
from sklearn.preprocessing import label_binarize
from sklearn import metrics

mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False

names = ['Age', 'Number of sexual partners', 'First sexual intercourse',
       'Num of pregnancies', 'Smokes', 'Smokes (years)',
       'Smokes (packs/year)', 'Hormonal Contraceptives',
       'Hormonal Contraceptives (years)', 'IUD', 'IUD (years)', 'STDs',
       'STDs (number)', 'STDs:condylomatosis',
       'STDs:cervical condylomatosis', 'STDs:vaginal condylomatosis',
       'STDs:vulvo-perineal condylomatosis', 'STDs:syphilis',
       'STDs:pelvic inflammatory disease', 'STDs:genital herpes',
       'STDs:molluscum contagiosum', 'STDs:AIDS', 'STDs:HIV',
       'STDs:Hepatitis B', 'STDs:HPV', 'STDs: Number of diagnosis',
       'STDs: Time since first diagnosis', 'STDs: Time since last diagnosis',
       'Dx:Cancer', 'Dx:CIN', 'Dx:HPV', 'Dx', 'Hinselmann', 'Schiller',
       'Citology', 'Biopsy']#df.columns
path = "datas/risk_factors_cervical_cancer.csv"  # 数据文件路径
data = pd.read_csv(path)

X = data[names[0:-4]]
Y = data[names[-4:]]
#随机森林可以处理多个目标变量的情况
X.head(5)
Y.head(5)

这个案例中需要预测的目标有四个目标值：Hiselmann、Schiller、Citlolgy、Biopsy。随机森林模型的一个特点是它可以同时预测多个属性。

下面的操作是对空值进行处理，之前讲过的案例中也有这步操作。但之前对空值的处理方式都是直接删除数据。但在实际工作中，样本数据可能不是很多。如果一发现异常就删除，可能导致样本量更少。

当特征是连续值的时候，可以考虑取对应特征的均值或去中位数填充。

如果特征是离散的，可以用众数来填充(谁多填谁)。

如果特征中有一个时间序列，那么缺失值对应的时间，和另一个时间如果相近，可以用另一个相近时间的样本对应的值填充。

__Imputer__函数，对应的单词是Data Imputation - 数据重组。这步操作涉及到统计学中的一个概念，当数据服从某种分布的时候，我们可以用某些策略对缺失的数据进行填充。

扩展：其中有一种叫MCMC的数据填充方法，有兴趣的读者可以自己查查。

#空值的处理
X = X.replace("?", np.NAN)
#使用Imputer给定缺省值，默认的是以mean
imputer = Imputer(missing_values="NaN")
X = imputer.fit_transform(X)

#数据分割
x_train,x_test,y_train,y_test = train_test_split(X, Y, 
    test_size=0.2, random_state=0)
print ("训练样本数量:%d,特征属性数目:%d,目标属性数目:%d" 
    % (x_train.shape[0],x_train.shape[1],y_train.shape[1]))
print ("测试样本数量:%d" % x_test.shape[0])

训练样本数量:686,特征属性数目:32,目标属性数目:4
测试样本数量:172

---

注意： 归一化本身不会影响数据的维度

#标准化
#分类模型，经常使用的是minmaxscaler归一化
#回归模型经常用standardscaler
ss = MinMaxScaler()
x_train = ss.fit_transform(x_train, y_train)
x_test = ss.transform(x_test)
x_train.shape

(686, 32)

---

PCA降维 - 文献：
https://www.cnblogs.com/pinard/p/6243025.html

#降维
# 降成1维
pca = PCA(n_components=1)
x_train = pca.fit_transform(x_train)
x_test = pca.transform(x_test)
x_train.shape
#解释了多少方差的比例
print(pca.explained_variance_ratio_)

[0.24021831]
解释了0.24，说明解释得不是很多。

---

n_estimators：建立100个模型
criterion：评价指标，这里用的是gini
max_depth：树的深度，这里设置1层，即森林分类一次结束。一般工作中设置5~10层。如果数据量非常大，可以再深一点。

#随机森林模型
forest = RandomForestClassifier(n_estimators=100, 
    criterion='gini', max_depth=1, random_state=0)
#max_depth一般不宜设置过大，把每个模型作为一个弱分类器
forest.fit(x_train, y_train)

RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',

        max_depth=1, max_features='auto', max_leaf_nodes=None,
        min_impurity_decrease=0.0, min_impurity_split=None,
        min_samples_leaf=1, min_samples_split=2,
        min_weight_fraction_leaf=0.0, n_estimators=100, n_jobs=1,
        oob_score=False, random_state=0, verbose=0, warm_start=False)

---

label_binarize：标签二值化
https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.label_binarize.html
形成哑编码：

标签二值化举例:
当我们的数据分类出现(0,1,2)个值，我们可以取0为正例，此时(1,2)为负例。由此将多分类转为二分类。

分析函数roc_curve的第一个参数:
label_binarize(y_test[names[-4]],classes=(0,1,2))[:,0:-1].ravel()

classes=(0,1,2) 代表有三个类别;
y_test[names[-4]] 表示Hislmann这一列; y_test 是真实值。
[:,0:-1] 除了最后一列之外都获取；为什么可以去掉一列？当k个不同数据形成k列哑编码时，其实可以用k-1列来代替。举个栗子：
A=100 B=010 C=001 如果去掉一列 A=10 B=01 C=00 同样可以即2(k-1)列哑编码可以表示3(k)个属性。

ravel() 用法举例：

总结:
label_binarize(y_test[names[-4]],classes=(0,1,2))[:,0:-1].ravel()

![label_binarize(y_test[names[-4]],classes=(0,1,2))[:,0:-1].ravel() - 1](https://upload-images.jianshu.io/upload_images/3153092-46a828269ca0c3ed.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

![label_binarize(y_test[names[-4]],classes=(0,1,2))[:,0:-1].ravel() - 2](https://upload-images.jianshu.io/upload_images/3153092-0832670fb9098bf3.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

![label_binarize(y_test[names[-4]],classes=(0,1,2))[:,0:-1].ravel() - 3](https://upload-images.jianshu.io/upload_images/3153092-bdf34857742384aa.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

---

分析函数roc_curve的第二个参数: forest_y_score[0].ravel() 预测值

forest_y_score = forest.predict_proba(x_test)# prodict_proba输出概率；

在计算roc_curve时，fpr和tpr分别是怎么算出来的？通过不同的阈值去判断概率。首先要

---

roc_curve：roc曲线评价指标
https://blog.csdn.net/Titan0427/article/details/79356290
https://blog.csdn.net/hh1294212648/article/details/77649127
该函数返回这三个变量：fpr,tpr,和阈值thresholds;
forest_fpr1 - fpr值、forest_tpr1 - tpr值、 _ - 这个值不想要，用下划线表示不接收。

ROC观察模型正确地识别正例的比例与模型错误地把负例数据识别成正例的比例之间的权衡。TPR的增加以FPR的增加为代价。ROC曲线下的面积是模型准确率的度量，AUC（Area under roccurve）。

__纵坐标：__真正率（True Positive Rate , TPR）或灵敏度（sensitivity）
TPR = TP /（TP + FN） （正样本预测结果数 / 正样本实际数）
__横坐标：__假正率（False Positive Rate , FPR）
FPR = FP /（FP + TN） （被预测为正的负样本结果数 /负样本实际数）

#模型效果评估
score = forest.score(x_test, y_test)
print ("准确率:%.2f%%" % (score * 100))
#模型预测
forest_y_score = forest.predict_proba(x_test)# prodict_proba输出概率
#计算ROC值
forest_fpr1, forest_tpr1, _ = metrics.roc_curve(
  label_binarize(y_test[names[-4]],classes=(0,1,2))[:,0:-1].ravel(), 
    forest_y_score[0].ravel())

forest_fpr2, forest_tpr2, _ = metrics.roc_curve(
    label_binarize(y_test[names[-3]],classes=(0,1,2))[:,0:-1].ravel(), 
      forest_y_score[1].ravel())

forest_fpr3, forest_tpr3, _ = metrics.roc_curve(
    label_binarize(y_test[names[-2]],classes=(0,1,2))[:,0:-1].ravel(), 
      forest_y_score[2].ravel())

forest_fpr4, forest_tpr4, _ = metrics.roc_curve(
    label_binarize(y_test[names[-1]],classes=(0,1,2))[:,0:-1].ravel(), 
        forest_y_score[3].ravel())
#AUC值
auc1 = metrics.auc(forest_fpr1, forest_tpr1)
auc2 = metrics.auc(forest_fpr2, forest_tpr2)
auc3 = metrics.auc(forest_fpr3, forest_tpr3)
auc4 = metrics.auc(forest_fpr4, forest_tpr4)

print ("Hinselmann目标属性AUC值：", auc1)
print ("Schiller目标属性AUC值：", auc2)
print ("Citology目标属性AUC值：", auc3)
print ("Biopsy目标属性AUC值：", auc4)

准确率:89.53%
Hinselmann目标属性AUC值： 0.9841806381828014
Schiller目标属性AUC值： 0.9497363439697133
Citology目标属性AUC值： 0.9376352082206598
Biopsy目标属性AUC值： 0.9520686857760952

---

label_binarize(y_test[names[-4]],classes=(0,1,2)).shape
# forest_y_score[0]
# label_binarize(y_test[names[-4]],classes=(0,1,2)).shape
# # y_test[names[-4]].value_counts()

(172, 3)

---

label_binarize(y_test[names[-4]],classes=(0,1,2))
y_test[names[-4]]
y_test.head(20)

---

##  画图（ROC图）
plt.figure(figsize=(8, 6), facecolor='w')
plt.plot(forest_fpr1,forest_tpr1,c='r',lw=2,label=u'Hinselmann目标属性,AUC=%.3f' % auc1)
plt.plot(forest_fpr2,forest_tpr2,c='b',lw=2,label=u'Schiller目标属性,AUC=%.3f' % auc2)
plt.plot(forest_fpr3,forest_tpr3,c='g',lw=2,label=u'Citology目标属性,AUC=%.3f' % auc3)
plt.plot(forest_fpr4,forest_tpr4,c='y',lw=2,label=u'Biopsy目标属性,AUC=%.3f' % auc4)
plt.plot((0,1),(0,1),c='#a0a0a0',lw=2,ls='--')
plt.xlim(-0.001, 1.001)
plt.ylim(-0.001, 1.001)
plt.xticks(np.arange(0, 1.1, 0.1))
plt.yticks(np.arange(0, 1.1, 0.1))
plt.xlabel('False Positive Rate(FPR)', fontsize=16)
plt.ylabel('True Positive Rate(TPR)', fontsize=16)
plt.grid(b=True, ls=':')
plt.legend(loc='lower right', fancybox=True, framealpha=0.8, fontsize=12)
plt.title(u'随机森林多目标属性分类ROC曲线', fontsize=18)
plt.show()

---

比较不同树数目、树最大深度的情况下随机森林的正确率。
一般情况下，初始的随机森林树个数是100，深度1，如果需要我们再进行优化操作。

x_train2,x_test2,y_train2,y_test2 = train_test_split(X, Y, test_size=0.5, random_state=0)
print ("训练样本数量%d，测试样本数量:%d" % (x_train2.shape[0], x_test2.shape[0]))
## 比较
estimators = [1,50,100,500]
depth = [1,2,3,7,15]
# x1, x2 = np.meshgrid(estimators, depth)
err_list = []
for es in estimators:
    es_list = []
    for d in depth:
        tf = RandomForestClassifier(n_estimators=es, criterion='gini', max_depth=d, max_features = None, random_state=0)
        tf.fit(x_train2, y_train2)
        st = tf.score(x_test2, y_test2)
        err = 1 - st
        es_list.append(err)
        print ("%d决策树数目，%d最大深度，正确率:%.2f%%" % (es, d, st * 100))
    err_list.append(es_list)

    
## 画图
plt.figure(facecolor='w')
i = 0
colors = ['r','b','g','y']
lw = [1,2,4,3]
max_err = 0
min_err = 100
for es,l in zip(estimators,err_list):
    plt.plot(depth, l, c=colors[i], lw=lw[i], label=u'树数目:%d' % es)
    max_err = max((max(l),max_err))
    min_err = min((min(l),min_err))
    i += 1
plt.xlabel(u'树深度', fontsize=16)
plt.ylabel(u'错误率', fontsize=16)
plt.legend(loc='upper left', fancybox=True, framealpha=0.8, fontsize=12)
plt.grid(True)
plt.xlim(min(depth),max(depth))
plt.ylim(min_err * 0.99, max_err * 1.01)
plt.title(u'随机森林中树数目、深度和错误率的关系图', fontsize=18)
plt.show()

训练样本数量429，测试样本数量:429
1决策树数目，1最大深度，正确率:86.48%
1决策树数目，2最大深度，正确率:86.95%
1决策树数目，3最大深度，正确率:84.62%
1决策树数目，7最大深度，正确率:82.75%
1决策树数目，15最大深度，正确率:78.09%
50决策树数目，1最大深度，正确率:86.71%
50决策树数目，2最大深度，正确率:86.48%
50决策树数目，3最大深度，正确率:86.48%
50决策树数目，7最大深度，正确率:86.25%
50决策树数目，15最大深度，正确率:84.38%
100决策树数目，1最大深度，正确率:86.95%
100决策树数目，2最大深度，正确率:86.25%
100决策树数目，3最大深度，正确率:86.48%
100决策树数目，7最大深度，正确率:86.25%
100决策树数目，15最大深度，正确率:85.08%
500决策树数目，1最大深度，正确率:86.48%
500决策树数目，2最大深度，正确率:86.48%
500决策树数目，3最大深度，正确率:86.48%
500决策树数目，7最大深度，正确率:86.25%
500决策树数目，15最大深度，正确率:84.85%