天天看點

07 內建學習 - 随機森林案例一:宮頸癌預測 - 預設資料填充政策、PCA降維、ROC曲線、标簽二值化

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() - 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%