GIS中如何实现随机森林算法源程序代码

一、问题描述

利用集成学习算法之一，随机森林实现宫颈癌数据集的多分类问题。

二、数据集分析

• 宫颈癌数据集下载地址为：http://archive.ics.uci.edu/ml/datasets/Cervical+cancer+(Risk+Factors)
• 宫颈癌数据集中一共有858个样本，共36列数据，目标有4个目标值：Hiselmann、Schiller、Citlolgy、Biopsy
• 数据集中有部分缺失值，用“？”标出
• 属性信息如下：
• (int) Age
  (int) Number of sexual partners
  (int) First sexual intercourse (age)
  (int) Num of pregnancies
  (bool) Smokes
  (bool) Smokes (years)
  (bool) Smokes (packs/year)
  (bool) Hormonal Contraceptives
  (int) Hormonal Contraceptives (years)
  (bool) IUD
  (int) IUD (years)
  (bool) STDs
  (int) STDs (number)
  (bool) STDs:condylomatosis
  (bool) STDs:cervical condylomatosis
  (bool) STDs:vaginal condylomatosis
  (bool) STDs:vulvo-perineal condylomatosis
  (bool) STDs:syphilis
  (bool) STDs:pelvic inflammatory disease
  (bool) STDs:genital herpes
  (bool) STDs:molluscum contagiosum
  (bool) STDs:AIDS
  (bool) STDs:HIV
  (bool) STDs:Hepatitis B
  (bool) STDs:HPV
  (int) STDs: Number of diagnosis
  (int) STDs: Time since first diagnosis
  (int) STDs: Time since last diagnosis
  (bool) Dx:Cancer
  (bool) Dx:CIN
  (bool) Dx:HPV
  (bool) Dx
  (bool) Hinselmann: target variable
  (bool) Schiller: target variable 席勒
  (bool) Cytology: target variable 细胞学
  (bool) Biopsy: target variable 活组织检查

三、代码实现

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl

from sklearn.preprocessing import Imputer
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import label_binarize
from sklearn.model_selection import train_test_split

from sklearn import tree
from sklearn.ensemble import RandomForestClassifier
from sklearn import metrics

import warnings
warnings.filterwarnings("ignore")

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']
data = pd.read_csv('risk_factors_cervical_cancer.csv')

# 查看数据
print(data.shape)

# 划分数据集中的特征和目标值
X = data[names[0:-4]]   # X = data.iloc[:,0:-4]
y = data[names[-4:]]    # Y = data.iloc[:,-4:]
# 预测的目标有四个目标值：Hiselmann、Schiller、Citlolgy、Biopsy
print(type(X))
print(X.head())
print(y.head())
# 空值处理
X = X.replace("?", np.NAN)
# 使用Imputer给定缺省值，默认以mean也就是均值替换
imputer = Imputer(missing_values="NaN")
X = imputer.fit_transform(X,y)

# 数据集分割
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])

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

# 模型构建
rf = RandomForestClassifier(n_estimators=50,criterion="gini",max_depth=1,\
                            random_state=10)
rf.fit(X_train,y_train)

# 模型效果评估
score = rf.score(X_test, y_test)
print("准确率：%.2f%%" %(score * 100))

#模型预测
rf_y_predict = rf.predict_proba(X_test)# prodict_proba输出概率
# print(rf_y_predict)

#计算ROC值
forest_fpr1, forest_tpr1, thresholds1 = \
             metrics.roc_curve(label_binarize(y_test[names[-4]],classes=(0,1,2)).\
                               T[0:-1].T.ravel(), rf_y_predict[0].ravel())
forest_fpr2, forest_tpr2, thresholds2 = \
             metrics.roc_curve(label_binarize(y_test[names[-3]],classes=(0,1,2)).\
                               T[0:-1].T.ravel(), rf_y_predict[1].ravel())
forest_fpr3, forest_tpr3, thresholds3 = \
             metrics.roc_curve(label_binarize(y_test[names[-2]],classes=(0,1,2)).\
                               T[0:-1].T.ravel(), rf_y_predict[2].ravel())
forest_fpr4, forest_tpr4, thresholds4 = \
             metrics.roc_curve(label_binarize(y_test[names[-1]],classes=(0,1,2)).\
                               T[0:-1].T.ravel(), rf_y_predict[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值：%.2f"%auc1)
print ("Schiller目标属性AUC值：%.2f"%auc2)
print ("Citology目标属性AUC值：%.2f"%auc3)
print ("Biopsy目标属性AUC值：%.2f"%auc4)

# 画图（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()

补充说明：

• Imputer类：该类在sklearn.preprocessing库中，用于处理丢失的数据，其格式为：

• imputer = Imputer(missing_values = "NaN", strategy = "mean", axis = 0)
• 参数意义：

• 1. missing_values：缺失值的占位符，所有出现的missing_values都将被估算，可以给它一个整数或者NaN来查找缺失值
• 2. strategy：插补策略，采用“mean”，则使用沿轴的平均值替换缺失值，还可以用“median”中位数和“most_frequent”最常见来替换
• 3. axis：可以指定0或者1，0沿列插入，1沿行插入

• imputer.fit()方法将数据相匹配，imputer.transform()方法将缺失值进行替换

为什么X_train标准化是用fit_transform()，而X_test标准化是用transform()呢？

• fit_transform()干了两件事：fit找到数据转换规则，并将数据标准化
• transform：是将数据进行转换，比如数据的归一化和标准化，将测试数据按照训练数据同样的模型进行转换，得到特征向量。可以直接把转换规则拿来用，所以并不需要fit_transform()，否则，两次标准化后的数据格式（或者说数据参数）就不一样了

• 归一化（Normalization）：属性缩放到一个指定的最大和最小值（通常是1-0）之间，可以通过preprocessing.MinMaxScaler类实现，常用的最小最大规范化方法(x-min(x))/(max(x)-min(x))
• 标准化（Standardization）：将数据按比例缩放，使之落入一个小的特定区间内，标准化后的数据可正可负，一般绝对值不会太大。计算时对每个属性/每列分别进行，将数据按期属性（按列进行）减去其均值，并处以其方差。得到的结果是，对于每个属性/每列来说所有数据都聚集在0附近，方差为1。使用z-score方法规范化(x-mean(x))/std(x)

四、随机森林算法

• 集成学习（Ensemble learning）通过构建并结合多个学习器来完成学习任务，有时也被成为多分类器系统、基于委员会的学习。一般结构是：先产生一组“个体学习器”，再用某种策略将它们结合起来，结合策略主要有平均法、投票法和学习法等。
• 要获得好的集成，个体学习器应“好而不同”，即个体学习器要有一定的“准确性”，即学习器不能太坏，并且要有多样性，即个体学习器间具有差异。
• 集成学习模型中所有的个体学习器如果都是一个种类的，叫做同质的，；如果不全是一种种类的，叫做异质的，同质个体学习器应用广泛，一般常说的集成学习方法都是指的同质个体学习器，而同质个体学习器使用最多的模型是CART决策树和神经网络。
• 同质个体学习器根据个体学习器的生成方式，将集成学习方法大致可以分为两类：

• 序列集成方法，其中参与训练的个体学习器按照顺序生成。个体学习器间存在强依赖关系，必须串行生成的序列化方法。通过对之前训练中错误标记的样本赋值较高的权重，可以提高整体的预测效果，代表为Boosting；
• 并行集成方法，其中参与训练的个体学习器并行生成，个体学习器间不存在强依赖关系、可同时生成的并行化方法，利用个体学习器之间的独立性，通过平均可以显著降低错误，代表为Bagging和随机森林

• 集成学习的组合策略：

• 平均法：对于数值类的回归预测问题，常用的组合策略是平均法，也就是对若干个弱学习器的输出进行平均得到最终的预测输出
• 投票法：对于分类问题的预测，常用的策略是投票法，最简单的投票法是相对多数投票法，也就是少数服从多数，稍微复杂的投票法是绝对多数投票法，也就是票数过半，更复杂的是加权投票法，每个弱学习器的分类票数乘以一个权重，最终将各个类别的加权票数求和，最大的值对应的类别为最终类别。
• 学习法：对于学习法，代表方法是stacking，当使用stacking的结合策略时， 我们不是对弱学习器的结果做简单的逻辑处理，而是再加上一层学习器，也就是说，我们将训练集弱学习器的学习结果作为输入，将训练集的输出作为输出，重新训练一个学习器来得到最终结果。

• 随机森林是集成学习的一种，是将多棵决策树整合成森林，并合起来预测最终结果。其中的每一棵树就是一个分类器，对于一个输入样本，N棵树就有N个分类结果，随机森林集成了所有的分类投票结果，将投票数最多的类别指定为最终结果，因此，这是一种Bagging的思想。
• Bagging的弱学习器之间没有直接联系，其特点在于”随机采样“（bootsrap），也就是从训练集里面采集固定个数的样本，每次采集一个样本后，都将其放回，之后采集样本的时候有可能再次被采集到。对于Bagging算法，一般会随机采集和训练集样本数m一样个数的样本。这样得到的采样集和训练集样本的个数相同，但是样本内容不同。
• 每棵树的生成规则：

• （1）如果训练集大小为N，对于每棵树而言，随机且有放回地从训练中抽取N个训练样本，作为该树的训练集，重复K次，生成K组训练样本集。
• （2）如果每个特征的样本维度为M，指定一个常数m<<M，随机地从M个特征中选取m个特征。
• （3） 利用m个特征对每棵树尽最大程度的生长，并且没有剪枝过程。

• 随机森林分类器：
• class sklearn.ensemble.RandomForestClassifier(n_estimators=10, criterion='gini', max_depth=None, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features=’auto’, max_leaf_nodes=None, min_impurity_decrease=0.0, min_impurity_split=None, bootstrap=True, oob_score=False, n_jobs=1, random_state=None, verbose=0, warm_start=False, class_weight=None)
• 主要参数：

• n_estimators : integer, optional (default=10)  整数，可选择(默认值为10)。
• criterion : string, optional (default=”gini”) 字符串，可选择(默认值为“gini”)。
• max_features : int, float, string or None, optional (default=”auto”)  整数，浮点数，字符串或者无值，可选的（默认值为"auto"），衡量分裂质量的性能（函数）。 受支持的标准是基尼不纯度的"gini",和信息增益的"entropy"（熵）。
• max_depth : integer or None, optional (default=None)  整数或者无值，可选的（默认为None），（决策）树的最大深度。如果值为None，那么会扩展节点，直到所有的叶子是纯净的，或者直到所有叶子包含少于min_sample_split的样本。