一、数据介绍
此数据是一份源于美国某经济学杂志上，分析研究波士顿房价( Boston House
Price)的数据集。数据集中的每一行数据都是对波士顿周边或城镇房价的描述：
CRIM: 城镇人均犯罪率
ZN: 住宅用地所占比例
INDUS: 城镇中非住宅用地所占比例
CHAS: CHAS 虚拟变量,用于回归分析
NOX: 环保指数
RM: 每栋住宅的房间数
AGE: 1940 年以前建成的自住单位的比例
DIS: 距离 5 个波士顿的就业中心的加权距离。
RAD: 距离高速公路的便利指数
TAX: 每一万美元的不动产税率
PRTATIO: 城镇中的教师学生比例
B: 城镇中的黑人比例
LSTAT: 地区中有多少房东属于低收入人群
MEDV: 自住房屋房价中位数（也就是均价）
二、任务介绍
1、通过数据挖掘对影响波士顿房价的因素进行分析。
2、搭建一个波士顿房价预测模型。
特征工程
数据预处理与特征分析
MEDV作为目标变量，也就是因变量，其他变量作为自变量。
import pandas as pdimport matplotlib.pyplot as pltimport matplotlibimport numpy as np
data = pd.read_csv('C:\\RAM\\housing.csv')
检查有没有数据中有没有空值
data.isnull().any().sum()
反馈是0，说明该数据不需要对其进行空值处理，可以放心的进行后续工作了。
查看各个特征的散点分布
pd.plotting.scatter_matrix(data, alpha=0.7, figsize=(10,10), diagonal='kde')
scatter_matrix
corr = data.corr()
结合相关性变量，大体了解下数据的整体分布。
特征选择
特征维度较大，为了保证模型的高效预测，需要进行特征选择。每种特征都有自己含义和数据量级，单纯地依靠方差来判断可能效果不好，直接使用与目标变量的相关性强的变量作为最终的特征变量。
通过相关系数法进行特征选择
x = data[['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT ']] y = data[['MEDV']]from sklearn.feature_selection import SelectKBestfrom sklearn.feature_selection import f_regression SelectKBest = SelectKBest(f_regression, k=3) bestFeature = SelectKBest.fit_transform(x,y) SelectKBest.get_support() x.columns[SelectKBest.get_support()]
这里我们看出和波士顿房价相关性最强的三个因素，分别是，RM(每栋住宅的房间数)，PTRATIO(城镇中的教师学生比例)，LSTAT(地区中有多少房东属于低收入人群)。
还是具备一定逻辑性的，首先，房子越大房价自然高（不管在哪个地域），其次，师生比与房价成反比，教育的重视，教育资源越是富裕的地方，生源就会大，师生比自然会降低，周边的房价会升高，这就是所谓的“学区房”概念吧。关于有多少房东属于低收入人群和房价的负相关关系，这个也比较好理解，各种原因吧。
features = data[['RM', 'PTRATIO', 'LSTAT ']] pd.plotting.scatter_matrix(features, alpha=0.7, figsize=(6,6), diagonal='hist')
特征归一化
from sklearn.preprocessing import MinMaxScaler
scaler = MinMaxScaler()for feature in features.columns:
features['标准化'+feature] = scaler.fit_transform(features[[feature]])#散点可视化，查看特征归一化后的数据font={ 'family':'SimHei'
}
matplotlib.rc('font', **font)
pd.plotting.scatter_matrix(features[['标准化RM', '标准化PTRATIO', '标准化LSTAT ']], alpha=0.7, figsize=(6,6), diagonal='hist')
模型选择与优化
数据集拆分
将数据拆分为训练数据及测试数据。
from sklearn.cross_validation import train_test_split x_train, x_test, y_train, y_test = train_test_split(features[['标准化RM', '标准化PTRATIO', '标准化LSTAT ']], y, test_size=0.3,random_state=33)
模型选择
采用交叉验证的方法对模型进行评估
from sklearn.model_selection import cross_val_predictfrom sklearn.model_selection import cross_val_score
房价预测方面，打算尝试以下方法：
线性回归
这应该是最简单也是最好理解的一种方法。使用支持向量回归模型SVR
之前仅仅是用SVM做过一些二分类问题（SVC），这次尝试下解决回归问题。KNN模型
思路上感觉KNN可以做，但是不确定效果，可以尝试下。决策树模型
和KNN一样，感觉是可以做，但是具体效果不太确定。
线性回归模型预测房价
from sklearn import linear_model lr = linear_model.LinearRegression() lr_predict = cross_val_predict(lr,x_train, y_train, cv=5) lr_score = cross_val_score(lr, x_train, y_train, cv=5) lr_meanscore = lr_score.mean()
SVR模型预测房价
尝试三种核，'linear', 'poly', 'rbf'
#SVRfrom sklearn.svm import SVR linear_svr = SVR(kernel = 'linear') linear_svr_predict = cross_val_predict(linear_svr, x_train, y_train, cv=5) linear_svr_score = cross_val_score(linear_svr, x_train, y_train, cv=5) linear_svr_meanscore = linear_svr_score.mean() poly_svr = SVR(kernel = 'poly') poly_svr_predict = cross_val_predict(poly_svr, x_train, y_train, cv=5) poly_svr_score = cross_val_score(poly_svr, x_train, y_train, cv=5) poly_svr_meanscore = poly_svr_score.mean() rbf_svr = SVR(kernel = 'rbf') rbf_svr_predict = cross_val_predict(rbf_svr, x_train, y_train, cv=5) rbf_svr_score = cross_val_score(rbf_svr, x_train, y_train, cv=5) rbf_svr_meanscore = rbf_svr_score.mean()
KNN模型
在KNN的回归模型当中，我们没法确定n_neighbors，因此需要最优化这个参数。
分别计算n_neighbors=[1,2,...,19,20]
from sklearn.neighbors import KNeighborsRegressor
score=[]for n_neighbors in range(1,21):
knn = KNeighborsRegressor(n_neighbors, weights = 'uniform' )
knn_predict = cross_val_predict(knn, x_train, y_train, cv=5)
knn_score = cross_val_score(knn, x_train, y_train, cv=5)
knn_meanscore = knn_score.mean()
score.append(knn_meanscore)
plt.plot(score)
plt.xlabel('n-neighbors')
plt.ylabel('mean-score')
image.png
从上图可以发现，随着n_neighbors的逐渐增大，模型预测能力逐渐增强，但是当n_neighbors大于2以后，模型评分趋于下降。因此我们选取n_neighbors=2。
n_neighbors=2 knn = KNeighborsRegressor(n_neighbors, weights = 'uniform' ) knn_predict = cross_val_predict(knn, x_train, y_train, cv=5) knn_score = cross_val_score(knn, x_train, y_train, cv=5) knn_meanscore = knn_score.mean()
决策树预测房价
和KNN类似，这里没法确定决策树的深度，因此令最大深度分别是1至10。
#Decision Treefrom sklearn.tree import DecisionTreeRegressor
score=[]for n in range(1,11):
dtr = DecisionTreeRegressor(max_depth = n)
dtr_predict = cross_val_predict(dtr, x_train, y_train, cv=5)
dtr_score = cross_val_score(dtr, x_train, y_train, cv=5)
dtr_meanscore = dtr_score.mean()
score.append(dtr_meanscore)
plt.plot(np.linspace(1,10,10), score)
plt.xlabel('max_depth')
plt.ylabel('mean-score')
发现当max_depth为[3, 4, 5]时，决策时模型评分处于极值的样子。
因此取max_depth为4。
n=4 dtr = DecisionTreeRegressor(max_depth = n) dtr_predict = cross_val_predict(dtr, x_train, y_train, cv=5) dtr_score = cross_val_score(dtr, x_train, y_train, cv=5) dtr_meanscore = dtr_score.mean()
模型评估
接下来，汇总下评分。
evaluating = { 'lr':lr_y_score, 'linear_svr':linear_svr_score, 'poly_svr':poly_svr_score, 'rbf_svr':rbf_svr_score, 'knn':knn_score, 'dtr':dtr_score
}
evaluating = pd.DataFrame(evaluating)
可视化
evaluating.plot.kde(alpha=0.6,figsize=(8,7))
所有模型的得分分布密度
evaluating.hist(color='k',alpha=0.6,figsize=(8,7))
所有模型的得分分布柱状图
从以上两张图发现，kernerl为poly的SVR得分受数据影响明显，而且得分偏低，其他的几个模型类似linear/rbf的SVR，dtr，都呈现出相同的趋势，KNN模型应该是算是截至在现在得分最高的模型。
模型优化
接下来想想办法，看看SVR还能不能被拯救。
首先对线性核进行最优化。
线性核我们通过更改惩罚系数C来查看对模型的影响。
lSVR_score=[]for i in [1,10,1e2,1e3,1e4]: linear_svr = SVR(kernel = 'linear', C=i) linear_svr_predict = cross_val_predict(linear_svr, x_train, y_train, cv=5) linear_svr_score = cross_val_score(linear_svr, x_train, y_train, cv=5) linear_svr_meanscore = linear_svr_score.mean() lSVR_score.append(linear_svr_meanscore) plt.plot(lSVR_score)
线性核随C的变化而变化
观察C为[1,10,100,1000]时对模型的影响，发现当C为10时，模型评分处于极值状态，随着C的增大，模型得分趋于稳定变化不大，因此认为C为10时模型最优。而sklearn关于线性核的默认状态C为1。
linear_svr = SVR(kernel = 'linear', C=10) linear_svr_predict = cross_val_predict(linear_svr, x_train, y_train, cv=5) linear_svr_score = cross_val_score(linear_svr, x_train, y_train, cv=5) linear_svr_meanscore = linear_svr_score.mean()
'poly'核优化，通过尝试修改惩罚系数C和degree。
polySVR_score=[]for i in [1,10,1e2,1e3,1e4]: for j in np.linspace(1,10,10): poly_svr = SVR(kernel = 'poly', C=i, degree=j) poly_svr_predict = cross_val_predict(poly_svr, x_train, y_train, cv=5) poly_svr_score = cross_val_score(poly_svr, x_train, y_train, cv=5) poly_svr_meanscore = poly_svr_score.mean() polySVR_score.append(poly_svr_meanscore) plt.plot(polySVR_score)
poly核最优化
从上图看出，对变量C来说，当C>10时，poly核模型得分普遍较高，当degree=2时，模型得分最高。我们查看SVR的API，发现，poly核默认的degree为3，C为1，可以从图中看出，优化后的模型得分增加不少。
poly_svr = SVR(kernel = 'poly', C=1000, degree=2) poly_svr_predict = cross_val_predict(poly_svr, x_train, y_train, cv=5) poly_svr_score = cross_val_score(poly_svr, x_train, y_train, cv=5) poly_svr_meanscore = poly_svr_score.mean()
‘rbf'核优化，优化惩罚系数C和gamma。
for i in [1,10,1e2,1e3,1e4]:
rbfSVR_score=[] for j in np.linspace(0.1,1,10):
rbf_svr = SVR(kernel = 'rbf', C=i, gamma=j)
rbf_svr_predict = cross_val_predict(rbf_svr, x_train, y_train, cv=5)
rbf_svr_score = cross_val_score(rbf_svr, x_train, y_train, cv=5)
rbf_svr_meanscore = rbf_svr_score.mean()
rbfSVR_score.append(rbf_svr_meanscore)
plt.plot(np.linspace(0.1,1,10),rbfSVR_score,label='C='+str(i))
plt.legend()
plt.xlabel('gamma')
plt.ylabel('score')
rbf核最优化
从上图发现，gamma从0.1递增到1.0，步长为0.1，模型得分有递增的趋势，当C>=10时，gamma对模型的影响较小，当C=100时，gamma=0.5时，’rbf‘模型评分更高。
rbf_svr = SVR(kernel = 'rbf', C=100, gamma=0.5) rbf_svr_predict = cross_val_predict(rbf_svr, x_train, y_train, cv=5) rbf_svr_score = cross_val_score(rbf_svr, x_train, y_train, cv=5) rbf_svr_meanscore = rbf_svr_score.mean()
接下来，对最优化的模型进行二次归类。
optimizer = { 'lr':lr_score, 'linear_svr':linear_svr_score, 'poly_svr':poly_svr_score, 'rbf_svr':rbf_svr_score, 'knn':knn_score, 'dtr':dtr_score
}
optimizer= pd.DataFrame(optimizer)可视化
optimizer.plot.kde(alpha=0.6,figsize=(8,7))
plt.xlabel('score')
最优化后的模型分数密度分布曲线
optimizer.hist(color='k',alpha=0.6,figsize=(8,7))
最优化后的模型得分分布柱状图
最优模型确定
对每个经过优化后的模型进行得分对比
此时发现，rbf核的SVR模型在优化后变成了最好的模型。线性核的SVR和线性回归因为策略的局限性，模型能力排在最后。
模型预测
接下来对测试数据集进行预测。
#rbfrbf_svr.fit(x_train,y_train)
rbf_svr_y_predict = rbf_svr.predict(x_test)
rbf_svr_y_predict_score=rbf_svr.score(x_test, y_test)#KNNknn.fit(x_train,y_train)
knn_y_predict = knn.predict(x_test)
knn_y_predict_score = knn.score(x_test, y_test)#poly_svrpoly_svr.fit(x_train,y_train)
poly_svr_y_predict = poly_svr.predict(x_test)
poly_svr_y_predict_score = poly_svr.score(x_test, y_test)#dtrdtr.fit(x_train, y_train)
dtr_y_predict = dtr.predict(x_test)
dtr_y_predict_score = dtr.score(x_test, y_test)#lrlr.fit(x_train, y_train)
lr_y_predict = lr.predict(x_test)
lr_y_predict_score = lr.score(x_test, y_test)#linear_svrlinear_svr.fit(x_train, y_train)
linear_svr_y_predict = linear_svr.predict(x_test)
linear_svr_y_predict_score = linear_svr.score(x_test, y_test)
predict_score = { 'lr':lr_y_predict_score, 'linear_svr':linear_svr_y_predict_score, 'poly_svr':poly_svr_y_predict_score, 'rbf_svr':rbf_svr_y_predict_score, 'knn':knn_y_predict_score, 'dtr':dtr_y_predict_score
}
predict_score = pd.DataFrame(predict_score, index=['score']).transpose()
predict_score.sort_values(by='score',ascending = False)预测结果排名
对各个模型的预测值整理
plt.scatter(np.linspace(0,151,152),y_test, label='predict data')
labelname=[ 'rbf_svr_y_predict', 'knn_y_predict', 'poly_svr_y_predict', 'dtr_y_predict', 'lr_y_predict', 'linear_svr_y_predict']
y_predict={ 'rbf_svr_y_predict':rbf_svr_y_predict, 'knn_y_predict':knn_y_predict[:,0], 'poly_svr_y_predict':poly_svr_y_predict, 'dtr_y_predict':dtr_y_predict, 'lr_y_predict':lr_y_predict[:,0], 'linear_svr_y_predict':linear_svr_y_predict
}
y_predict=pd.DataFrame(y_predict)for name in labelname:
plt.plot(y_predict[name],label=name)
plt.xlabel('predict data index')
plt.ylabel('target')
plt.legend()
综合表现来看，rbf_svr的鲁棒性更强一点。
2018/09/01
关于特征选择过程的思考：在上述步骤中，通过将特征与目标变量进行相关系数并作为依据，选出三个具有强相关性的特征作为模型的输入。但是似乎忽略了弱相关变量对目标变量的影响。将在以下部分进行讨论。
各个特征之间虽然没有很强的一对一字面上理解的因果关系，但是，根据各个单一特征之间的皮尔森相关系数，尝试排除彼此有很强线性相关的变量，保证特征之间都是相互独立非线性因果的关系。这样做的目的也是为了避免多重共线性的过拟合问题。
for column in corr.columns: inxs = (abs(corr[column])<1) & (abs(corr[column])>0.6) for inx in corr.index[inxs]: print(column+'<-->'+inx) corr['MEDV'].abs().sort_values(ascending=False)
这时候输出存在线性强相关的逻辑变量。
CRIM<-->TAX
CRIM<-->NOX
CRIM<-->RAD
ZN<-->DIS
INDUS<-->LSTAT
INDUS<-->TAX
INDUS<-->NOX
INDUS<-->AGE
INDUS<-->DIS
NOX<-->INDUS
NOX<-->TAX
NOX<-->RAD
NOX<-->AGE
NOX<-->CRIM
NOX<-->DIS
RM<-->LSTAT
RM<-->MEDV
AGE<-->LSTAT
AGE<-->INDUS
AGE<-->NOX
AGE<-->DIS
DIS<-->INDUS
DIS<-->NOX
DIS<-->AGE
DIS<-->ZN
RAD<-->TAX
RAD<-->NOX
RAD<-->CRIM
TAX<-->INDUS
TAX<-->NOX
TAX<-->RAD
TAX<-->CRIM
LSTAT <-->INDUS
LSTAT <-->AGE
LSTAT <-->RM
LSTAT <-->MEDV
MEDV<-->LSTAT
MEDV<-->RM
根据各个特征与目标变量的相关系数绝对值的排名，依次对强相关的特征进行最优筛选。
最终获得的特征为：
['LSTAT','PTRATIO','TAX','ZN','B','CHAS']
以上六个特征是与目标变量最相关最独立的几个特征。
总的特征我们在以上基础上，将其余的按照与目标变量的相关系数依次导入。
['LSTAT','PTRATIO','TAX','ZN','B','CHAS','RM','INDUS','NOX','RAD','AGE','CRIM','DIS']
在这里，将最少特征数量设定为1，最大特征数设为已知的6，其实这也算是一种贪心算法，仅作为探索的尝试吧。
通过上边的结果，我们选取效果比较好的几个模型，决策树模型，knn模型，和'rbf'核的SVR模型。
knn_meanscore = []
dtr_meanscore = []
rbf_svr_meanscore = []new = features[['标准化LSTAT ','标准化PTRATIO','标准化TAX','标准化ZN','标准化B','标准化CHAS','标准化RM','标准化INDUS','标准化NOX','标准化RAD','标准化AGE','标准化CRIM','标准化DIS']]for m in range(1,14): print('feature num: %d\r\n' %m)
selectKBest = SelectKBest(f_regression, k=int(m))
bestFeature = selectKBest.fit_transform(new,y)
xx = new[new.columns[selectKBest.get_support()]] print("features' name: ", xx.columns.tolist())
x_train, x_test, y_train, y_test = train_test_split(xx,y,test_size=0.3,random_state=33) #knn
n_neighbors=2
knn = KNeighborsRegressor(n_neighbors, weights = 'uniform' )
knn_score = cross_val_score(knn, x_train, y_train, cv=5)
knn_meanscore.append(knn_score.mean()) #Decision Tree
dtr = DecisionTreeRegressor(max_depth = 4)
dtr_score = cross_val_score(dtr, x_train, y_train, cv=5)
dtr_meanscore.append(dtr_score.mean()) #rbf
rbf_svr = SVR(kernel = 'rbf',C=100,gamma=0.5)
rbf_svr_score = cross_val_score(rbf_svr, x_train, y_train, cv=5)
rbf_svr_meanscore.append(rbf_svr_score.mean())
evaluating = { 'rbf_svr':rbf_svr_meanscore, 'knn':knn_meanscore, 'dtr':dtr_meanscore
}
evaluating = pd.DataFrame(evaluating)
evaluating.plot()
plt.xticks(evaluating.index,range(1, 13))
plt.xlabel("features' number")
plt.ylabel('Model score')
从以上可以看出三种模型在随着排序好的特征数逐渐增加的过程中，模型得分逐渐趋于稳定。
实际相互独立的特征只有6个，我们看出，从第7个特征开始，模型的评分就趋于稳定，这和我们开始假设是吻合的，冗余的特征对于模型并没有特别的贡献。对于决策树来说，似乎特征数量愈多，模型得分越高。
对于rbf核的SVR来说，它的得分是最高的，0.83，这比之前只选取强相关的三个特征的模型得分高出10%。
接下来，我们确定特征数为6，对test数据进行预测。
selectKBest = SelectKBest(f_regression, k=6)
bestFeature = selectKBest.fit_transform(new,y)
xx = new[new.columns[selectKBest.get_support()]]print("features' name: ", xx.columns.tolist())
#knnn_neighbors=2knn = KNeighborsRegressor(n_neighbors, weights = 'uniform' )
knn.fit(x_train, y_train)
knn_y_predict = knn.predict(x_test)
knn_y_predict_score = knn.score(x_test, y_test)#Decision Treedtr = DecisionTreeRegressor(max_depth = 4)
dtr.fit(x_train, y_train)
dtr_y_predict = dtr.predict(x_test)
dtr_y_predict_score = dtr.score(x_test, y_test)#rbfrbf_svr = SVR(kernel = 'rbf',C=100,gamma=0.5)
rbf_svr.fit(x_train,y_train)
rbf_svr_y_predict = rbf_svr.predict(x_test)
rbf_svr_y_predict_score=rbf_svr.score(x_test, y_test)
predict_score = { 'rbf_svr':rbf_svr_y_predict_score, 'knn':knn_y_predict_score, 'dtr':dtr_y_predict_score
}
predict_score = pd.DataFrame(predict_score, index=['score']).transpose()
predict_score.sort_values(by='score',ascending = False)
可以看出rbf核的SVR比之前的模型预测得分高了3%。
其实运算的时间成本上也需要考虑，这里先不写了。
作者：要转算法的巨巨
链接：https://www.jianshu.com/p/9d3c51cf32b8
随时随地看视频
