一、数据介绍
此数据是一份源于美国某经济学杂志上,分析研究波士顿房价( 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 pd
import matplotlib.pyplot as plt
import matplotlib
import 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 SelectKBest
from 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_predict
from 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'
#SVR
from 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')
从上图可以发现,随着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 Tree
from 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为[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)
从上图看出,对变量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')
从上图发现,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和线性回归因为策略的局限性,模型能力排在最后。
-
模型预测
接下来对测试数据集进行预测。
#rbf
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)
#KNN
knn.fit(x_train,y_train)
knn_y_predict = knn.predict(x_test)
knn_y_predict_score = knn.score(x_test, y_test)
#poly_svr
poly_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)
#dtr
dtr.fit(x_train, y_train)
dtr_y_predict = dtr.predict(x_test)
dtr_y_predict_score = dtr.score(x_test, y_test)
#lr
lr.fit(x_train, y_train)
lr_y_predict = lr.predict(x_test)
lr_y_predict_score = lr.score(x_test, y_test)
#linear_svr
linear_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())
#knn
n_neighbors=2
knn = 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 Tree
dtr = 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)
#rbf
rbf_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%。
其实运算的时间成本上也需要考虑,这里先不写了。