数据挖掘入门指南：以kaggle:titanic为例

一般的数据挖掘竞赛或者项目包括以下步骤：

（1）数据预处理和特征工程

（2）构造模型

（3）模型融合

这篇博客将以kaggle上的titanic题目作为例子，具体讲一讲各个步骤的操作过程。数据集可自行到
kaggle:titanic 下 载。

PS：对于titanic这道题，个人觉得重点不是在于把成绩提高到多少，因为样本实在太少了，所以我们应该把关注点放在学习如何分析数据，以及构建模型等。

（1）数据预处理和特征工程

一般情况下，我们拿到的数据不可能是完全符合要求的，存在一些缺失值或者离群点之类的情况。

先读入数据：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
#读入数据
train = pd.read_csv('./train.csv')
test = pd.read_csv('./test.csv')
survived = train['Survived']
combine = pd.concat([train.drop('Survived',axis=1),test],axis=0) ##合并训练集和测试集

看一下缺失值情况：

combine.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 1309 entries, 0 to 417
Data columns (total 11 columns):
PassengerId    1309 non-null int64
Pclass         1309 non-null int64
Name           1309 non-null object
Sex            1309 non-null object
Age            1046 non-null float64
SibSp          1309 non-null int64
Parch          1309 non-null int64
Ticket         1309 non-null object
Fare           1308 non-null float64
Cabin          295 non-null object
Embarked       1307 non-null object
dtypes: float64(2), int64(4), object(5)
memory usage: 122.7+ KB

#也可以使用combine.isnull().sum()一目了然
combine.isnull().sum()
Out[356]:
PassengerId       0
Pclass            0
Name              0
Sex               0
Age             263
SibSp             0
Parch             0
Ticket            0
Fare              1
Cabin          1014
Embarked          2
dtype: int64

可以看到，Age,Fare,Cabin,Embarked是存在缺失值的，需要我们在后续补全缺失值。

接下来，使用matplotlib库和seaborn库画图对特征进行分析：

首先，将用户分成存活乘客和非存活乘客：

surv = train[train['Survived']==1]
nosurv = train[train['Survived']==0]

画连续特征可以用sns.distplot()函数，画离散特征可以用sns.barplot()函数：

surv_color = 'blue'
nosurv_color = 'red'

##Age:画年龄分布柱状图
plt.figure(figsize=[12,10])
plt.subplot(331)
sns.distplot(surv['Age'].dropna().values,bins=range(0,80,1),kde=False,color=surv_color,axlabel='surv_Age')
plt.subplot(332)
sns.distplot(nosurv['Age'].dropna().values,bins=range(0,80,1),kde=False,color=nosurv_color,axlabel='nosurv_Age')
##在某个特征下的幸存率
##eg:不同性别下的幸存率分布
plt.subplot(333)
sns.barplot('Sex','Survived',data=train)
plt.subplot(334)
sns.barplot('Pclass','Survived',data=train)
plt.subplot(335)
sns.barplot('Embarked','Survived',data=train)
plt.subplot(336)
sns.barplot('Parch','Survived',data=train)
plt.subplot(337)
sns.barplot('SibSp','Survived',data=train)
##因为Fare是连续值且分布范围很广，可以先对Fare取对数后再进行显示
plt.subplot(338)
sns.distplot(np.log10(surv['Fare'].dropna().values+1),kde=False,color=surv_color)
sns.distplot(np.log10(nosurv['Fare'].dropna().values+1),kde=False,color=nosurv_color,axlabel='Fare')
##统计乘客的家庭成员个数并显示不同家庭成员数量下的存活率
train['Family'] = train['SibSp'] + train['Parch']
plt.subplot(339)
sns.barplot('Family','Survived',data=train)

由图Family-Survived（即最后一个图）我们可以知道家庭成员个数在1，2,3的乘客存活率较高，所以我们可以把它区分开来：

combine['Family'] = combine['SibSp'] + combine['Parch']
combine['FamilyBins'] = np.where(combine['Family']==0,1,np.where(combine['Family']<4,0,1)) ##大于0小于4的值设为0，其他设为1

我们由前面知道特征'Age'存在很多缺失值，使用corr()看一下与年龄相关的特征有哪些？使用sns.heatmap()函数画图：

##显示各特征之间的相关性
plt.figure(figsize=(12,10))
corr = sns.heatmap(train.drop('PassengerId',axis=1).corr(),vmax=0.8,annot=True)

可以看到，与年龄相关性比较大的特征有：‘Pclass’，‘SibSp’，'Parch'

接下来，补全缺失值：

group = combine.groupby(['Pclass','SibSp','Parch']).Age
combine['Age'] = group.transform(lambda x:x.fillna(x.median())) #transform会将一个函数应用到各个分组，然后将结果放在适当的位置.

探索在不同性别下，用户在不同年龄下的生存情况：

##分离出不同性别生存与否的乘客
msurv = train[(train['Survived']==1) & (train['Sex']=='male')]
mnosurv = train[(train['Survived']==0) & (train['Sex']=='male')]
fsurv = train[(train['Survived']==1) & (train['Sex']=='female')]
fnosurv = train[(train['Survived']==0) & (train['Sex']=='female')]

分别显示在不同性别下，乘客在不同年龄下的幸存率：

plt.figure(figsize=[12,10])
plt.subplot(221)
sns.distplot(msurv['Age'].dropna().values,bins=range(0,80,1),kde=False,color=surv_color)
sns.distplot(mnosurv['Age'].dropna().values,bins=range(0,80,1),kde=False,color=nosurv_color,axlabel='maleAge')
plt.subplot(222)
sns.distplot(fsurv['Age'].dropna().values,bins=range(0,80,1),kde=False,color=surv_color)
sns.distplot(fnosurv['Age'].dropna().values,bins=range(0,80,1),kde=False,color=nosurv_color,axlabel='femaleAge')

由上图可以看到，在青中年阶段（大概在18-40岁），男性乘客的存活率明显偏低，女性乘客的存活率明显偏高。

因此，可以构造新特征：将年龄在18-40岁的男乘客归为0，女乘客归为1，其他归为2：

age_male_name = combine[((combine['Age']>=18) & (combine['Age']<40)) & (combine['Sex']=='male')]['Name'].values
age_female_name = combine[((combine['Age']>=18) & (combine['Age']<40)) & (combine['Sex']=='female')]['Name'].values
combine['AgeClass'] = np.where(combine['Name'].isin(age_male_name),0,np.where(combine['Name'].isin(age_female_name),1,2))

如果，我们要统计不同Pclass，各年龄情况的幸存率（年龄为近似连续值），该怎么画图呢？

可以使用sns.violinplot()函数：

sns.violinplot(x='Pclass',y='Age',hue='Survived',data=train,split=True)  #split=True能够把每个Pclass下，取Survived=0和Survived=1每个图的左右各一半
#合在一起，方便观察；
plt.hlines([0,12],xmin=-1,xmax=3,linestyles='dotted') #这一句是标定两条分界线（图中虚线），方便观察

PS：有一个小问题，应该是图片被稍微拉伸了，年龄是不会有负值的；但这不影响我们对整个分布的观察。

由上图可以看到，在年龄比较小时（图中小于12），Pclass=1,2的孩子基本上都存活下来了，Pclass=3时有部分孩子没有存活。

统计不同港口（Embarked），不同Pclass下男性和女性的幸存率：

sns.factorplot(x='Pclass',y='Survived',hue='Sex',col='Embarked',data=train)

如果想将折线图转换为柱状图，可以加入kind参数：kind='bar':

sns.factorplot(x='Pclass',y='Survived',hue='Sex',col='Embarked',data=train,kind='bar')

我们可以看到，Pclass=3时，Sex=male的乘客没有存活的，Sex=female的乘客全部存活，有异常吗？我们分析一下：

train[(train['Pclass']==1) & (train['Embarked']=='Q')][['Sex','Survived']]
Out[441]:
Sex  Survived
245    male         0
412  female         1

train[(train['Pclass']==2) & (train['Embarked']=='Q')][['Sex','Survived']]
Out[442]:
Sex  Survived
303  female         1
322  female         1
626    male         0

通过分析，我们可以得到，出现上述情况的原因是：在Embarked=Q上船且Pclass=1,2的乘客非常少，所以出现上图情况很正常。

构造新特征：将上图上幸存率较高和较低的情况分开，构成新的特征：

PSM_name = combine[((combine['Pclass']<3) & (combine['Sex']=='female')) | ((combine['Pclass']==3) & (combine['Sex']=='female') & (combine['Embarked']!='S'))]['Name'].values
combine['PSM'] = np.where(combine['Name'].isin(PSM_name),0,1)   #因为乘客中没有出现重名的情况，这里借用了名字的唯一性，使用Name作为中间变量

统计不同港口下，各Pclass等级乘客的比例，使用pd.crosstab()构造表格：

tab = pd.crosstab(combine['Embarked'],combine['Pclass'])
pic = tab.div(tab.sum(1).astype(float),axis=0).plot(kind='bar',stacked=True)
pic = plt.xlabel('Embarked')
pic = plt.ylabel('Percent')

tab
Out[465]:
Pclass      1    2    3
Embarked
C         141   28  101
Q           3    7  113
S         177  242  495

接下来，讲一下Fare特征要怎么画图分析。

因为Fare的数值范围比较大[0:600]，所以我们可以先给Fare取对数，再进行分析：

plt.figure(figsize=[12,10])
plt.subplot(311)
sns.distplot(np.log10(surv['Fare'].dropna().values+1),kde=False,color=surv_color)
sns.distplot(np.log10(nosurv['Fare'].dropna().values+1),kde=False,color=nosurv_color,axlabel='Fare')

可以看到，Fare越大，乘客幸存率越高。

统计不同Pclass等级下个Fare水平的幸存率，使用sns.boxplot()函数：

sns.boxplot(x='Pclass',y='Fare',hue='Survived',data=train).set_yscale('log')

补全缺失值：

1. 补全Fare的缺失值：由上述相关关系表知道与Fare相关性较大的特征有Pclass,Parch,SibSp，因此：

##fill nan values in Fare
nullFares = combine[combine.Fare.isnull()].index.values

combine.loc[nullFares,'Fare'] = combine[(combine['Pclass'] == 3) & (combine['Parch'] == 0) & (combine['SibSp'] == 0)].Fare.median()

2.补全Embarked的缺失值：

##查明缺失值的信息
combine[combine['Embarked'].isnull()]
Out[578]:
PassengerId  Pclass                                       Name     Sex  \
61            62       1                        Icard, Miss. Amelie  female
829          830       1  Stone, Mrs. George Nelson (Martha Evelyn)  female

Age  SibSp  Parch  Ticket  Fare Cabin Embarked
61   38.0      0      0  113572  80.0   B28      NaN
829  62.0      0      0  113572  80.0   B28      NaN

##与缺失值相同特征的乘客数量
combine.where((combine['Pclass']==1) & (combine['Sex']=='female')).groupby(['Embarked','Pclass','Sex','Parch','SibSp']).size()
Out[577]:
Embarked  Pclass  Sex     Parch  SibSp
C         1.0     female  0.0    0.0      30
1.0      20
1.0    0.0      10
1.0       6
2.0    0.0       2
2.0       2
3.0    1.0       1
Q         1.0     female  0.0    1.0       2
S         1.0     female  0.0    0.0      20
1.0      20
2.0       3
1.0    0.0       7
1.0       6
2.0    0.0       4
1.0       5
3.0       3
4.0    1.0       1

由上述统计可知，可使用Embarked='C'补全缺失值：

##fill nan values in Embarked
nullEmbarkeds = combine[combine.Embarked.isnull()].index.values
combine['Embarked'].iloc[nullEmbarkeds] = 'C'

（2）构建模型

这里使用xgboost的分类算法

from xgboost.sklearn import XGBClassifier
from sklearn.model_selection import cross_val_score
########################### offline model ###########################
##交叉验证得到线下训练的准确率
model = XGBClassifier(max_depth=6, n_estimators=1000, learning_rate=0.01)
scores = cross_val_score(model,x_train,y_train,cv=3)
print ('accuracy:{0:.5f}'.format(np.mean(scores)))

##使用xgboost的get_fscore得到特征的重要性并排序
model.fit(x_train, y_train)
importance = model.booster().get_fscore()
sort_importance = sorted(importance.items(),key=operator.itemgetter(1),reverse=True)
df = pd.DataFrame(sort_importance, columns=['feature', 'fscore'])
df['fscore'] = df['fscore'] / df['fscore'].sum()
print (df)

结果显示如下：



（3）模型融合

这里主要讲两层的模型融合，先简单说明一下思路（以二折交叉验证为例）：

1.将训练样本平均分为两部分：A1，A2；测试样本记为C；

2.用A1训练，A2测试，得到测试结果B2; 用A1训练，C测试，得到测试结果D2;

3.用A2训练，A1测试，得到测试结果B1; 用A2训练，C测试，得到测试结果D1;

4.由B1，B2组成第二层的训练集，由D1，D2组成第二层的测试集，标签为样本原来的真实特征；

5.使用第4步得到的数据集继续做预测（可以用逻辑回归等模型），得到预测结果。

以Titanic为例，代码如下：

from sklearn.ensemble import RandomForestClassifier,AdaBoostClassifier,GradientBoostingClassifier,ExtraTreesClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.cross_validation import KFold
from xgboost.sklearn import XGBClassifier
from sklearn.svm import SVC

NKFOLD = 5
len_train = train.shape[0]
len_test = test.shape[0]
kf = KFold(len_train,n_folds=5,random_state=0)

class Classifier(object):
def __init__(self,clf,param=None):
self.clf = clf(**param)

def train(self,x_train,y_train):
self.clf.fit(x_train,y_train)

def predicted(self,x_test):
return self.clf.predict(x_test)

##构建第一层模型
def getFirstModel(clf,x_train,y_train,x_test):
train_pred = np.zeros((len_train,))
test_pred = np.zeros((NKFOLD,len_test))
test_pred_mean = np.zeros((len_test,))
for i,(train_index,test_index) in enumerate(kf):
x_tr = x_train[train_index]
y_tr = y_train[train_index]
x_te = x_train[test_index]

clf.train(x_tr,y_tr)
train_pred[test_index] = clf.predicted(x_te)
test_pred[i,:] = clf.predicted(x_test)
test_pred_mean[:] = np.mean(test_pred,axis=0)
return train_pred.reshape((-1,1)),test_pred_mean.reshape((-1,1))

##各模型参数
rf_params = {
'n_jobs': -1,
'n_estimators': 500,
'warm_start': True,
#'max_features': 0.2,
'max_depth': 6,
'min_samples_leaf': 2,
'max_features' : 'sqrt',
'verbose': 0,
'random_state':0
}

# Extra Trees Parameters
et_params = {
'n_jobs': -1,
'n_estimators':500,
#'max_features': 0.5,
'max_depth': 8,
'min_samples_leaf': 2,
'verbose': 0,
'random_state':0
}

# AdaBoost parameters
ada_params = {
'n_estimators': 500,
'learning_rate' : 0.75,
'random_state':0
}

# Gradient Boosting parameters
gb_params = {
'n_estimators': 500,
#'max_features': 0.2,
'max_depth': 5,
'min_samples_leaf': 2,
'verbose': 0,
'random_state':0
}

xgb_params = {
'n_estimators' : 2000,
'max_depth': 4,
'min_child_weight': 2,
#gamma=1,
'gamma':0.9,
'subsample':0.8,
'colsample_bytree':0.8,
'objective': 'binary:logistic',
'nthread': -1,
'scale_pos_weight':1,
'seed':0
}

# Support Vector Classifier parameters
svc_params = {
'kernel' : 'linear',
'C' : 0.025,
'random_state':0
}

rf = Classifier(RandomForestClassifier,param = rf_params)
ext = Classifier(ExtraTreesClassifier,param = et_params)
ada = Classifier(AdaBoostClassifier,param = ada_params)
gbdt = Classifier(GradientBoostingClassifier,param = gb_params)
xgb = Classifier(XGBClassifier,param = xgb_params)
svc = Classifier(SVC,param = svc_params)

rf_train_pred,rf_test_pred = getFirstModel(rf,x_train,y_train,x_test)
ext_train_pred,ext_test_pred = getFirstModel(ext,x_train,y_train,x_test)
ada_train_pred,ada_test_pred = getFirstModel(ada,x_train,y_train,x_test)
gbdt_train_pred,gbdt_test_pred = getFirstModel(gbdt,x_train,y_train,x_test)
xgb_train_pred,xgb_test_pred = getFirstModel(xgb,x_train,y_train,x_test)
svc_train_pred,svc_test_pred = getFirstModel(svc,x_train,y_train,x_test)

x_train = np.concatenate((rf_train_pred,ext_train_pred,ada_train_pred,gbdt_train_pred,xgb_train_pred,svc_train_pred),axis=1)
x_test = np.concatenate((rf_test_pred,ext_test_pred,ada_test_pred,gbdt_test_pred,xgb_test_pred,svc_test_pred),axis=1)

model = LogisticRegression().fit(x_train,y_train)
prediction = model.predict(x_test)
result = pd.DataFrame({'PassengerId': PassengerId,
'Survived': prediction
})
result.to_csv('./result.csv',index=False)

补充一点：

在模型融合中，判断第一层某一个模型的预测值是否要加入到stacking中（第二层模型以线性回归为例）：

（1）将所有模型的预测值作为输入，由线性回归计算验证集的误差；

（2）对于第一层中的每一个模型，计算移除该模型后，混叠模型由线性回归计算出来的误差；

（3）移除对误差缩小贡献最小的一个模型，并回到第一步继续工作。

典型的，我们把贡献值低于(2-3)*10^-6的模型排除掉（值大小具体情况具体分析）

参考：

pytanic  （该kernel的作者写的非常不错，我是按照他的逻辑整理的）

模型融合参考自：Introduction to Ensembling/Stacking in Python

kaggle上关于模型融合的介绍：Kaggle Ensembling Guide