kaggle之房价预测

题目简述

House Prices: Advanced Regression Techniques
给出了可能和房价有关的79个变量,使用回归的方法预测房价
主要按照以下文章的思路,注意这篇文章只是做一个简单的笔记整理,供后查阅,因此是针对博主而言的:
数据分析:Comprehensive data exploration with Python
数据处理及线性回归在训练集上的比较:A study on Regression applied to the Ames dataset
xgboost+线性回归:Regularized Linear Models
stacking的使用:Stacked Regressions : Top 4% on LeaderBoard
最后简单的处理了一下数据成绩为0.11510
代码放在github上:

数据分析

这篇文章对数据的分析可谓经典,从中学到很多,但也有很多不是太明白
使用describe来判断DataFrame列的个数,均值,标准差,最小值,25%分位点,75%分位点,最大值

1
df_train['SalePrice'].describe()

1
2
3
4
5
6
7
8
9
count 1460.000000
mean 180921.195890
std 79442.502883
min 34900.000000
25% 129975.000000
50% 163000.000000
75% 214000.000000
max 755000.000000
Name: SalePrice, dtype: float64

关于数据的偏度和峰度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#偏度,峰度
#参见:
# http://blog.csdn.net/lxb1022/article/details/76919362
# http://blog.sina.com.cn/s/blog_514922200101c2dn.html
"""
偏度是描述数据分布形态的统计量,其描述的是某总体取值分布的对称性,
偏度大于0表示其数据分布形态与正太分布相比为正偏或右偏;有一条长尾巴在右边;反之左偏,即有一条长尾巴在左边
峰度是描述总体中所有取值分布形态陡缓程度的统计量。
这个统计量需要与正态分布相比较,峰度为0表示该总体数据分布与正态分布的陡缓程度相同;峰度大于0表示该总体数据分布与正态分布相比较为陡峭,为尖顶峰;
峰度小于0表示该总体数据分布与正态分布相比较为平坦,为平顶峰。峰度的绝对值数值越大表示其分布形态的陡缓程度与正态分布的差异程度越大。
"""
#价格变量的偏差和峰度,
print("Skewness: %f" % df_train['SalePrice'].skew())
print("Kurtosis: %f" % df_train['SalePrice'].kurt())

散点图,绘制两个变量之间的关系,concat为连接,axis=1表示列上的连接

1
2
3
4
#scatter画散点图,下面是GrLivArea和SalePrice之间的关系图
var='GrLivArea'
data=pd.concat([df_train['SalePrice'],df_train[var]],axis=1)#axis表示在列上进行连接
data.plot.scatter(x=var,y='SalePrice',ylim=(0,800000));

kaggle0101

箱形图的绘制,由异常点,分位点组成,可以观察数据的集中程度,异常情况等,Plotting with categorical data中给出的解释是说散点图可以很好的绘制出两个连续变量(数值变量)之间的关系,但是对于类别变量和数值变量之间的关系就需要用到类似的箱型图了,无意中发现了知乎专栏对seaborn全部的翻译Seaborn(sns)官方文档学习笔记
下图中,我就不懂这个异常点是怎么回事,如果按照给出的文章箱形图中异常点的检测来看异常点的话,计算公式为Q1-k(Q3-Q1)和Q3+k(Q3-Q1),那么异常点应该离Q1点比较远的距离吧,但是图上为什么不是这样?

1
2
3
4
5
var='OverallQual'
data=pd.concat([df_train['SalePrice'],df_train[var]],axis=1)
f,ax=plt.subplots(figsize=(8,6))#创建800*600像素的图形,返回figure图像和一个子图ax的列表
fig=sns.boxplot(x=var,y='SalePrice',data=data)
fig.axis(ymin=0,ymax=800000);#设置y的范围

kaggle0102
DataFrame中生成相关系数矩阵,sns通过热力图展现

1
2
3
4
5
#将生成的相关系数矩阵以热力图的形式展现出来
#这个corr是计算DataFrame中列之间的相关系数,如果存在非数字列,则会自动排除
corrmat = df_train.corr()#相关系数,比如有n列,则相关系数矩阵大小为n*n,其中的(i,j)表示的是第i列和第j列之间的相关系数,最后生成的是DataFrame
f, ax = plt.subplots(figsize=(12, 9))
sns.heatmap(corrmat, vmax=.8, square=True);#vmax=.8表示设置最大值为0.8;square表示将Axes横纵比设置为相等,这样每个单元格都是正方形

kaggle0103
上图可以清楚的显示有两个明显的淡黄块,即1stFlrSF与TotalBsmtSF,GarageCars和GarageArea,它们之间的相关性最强
这里提一下方差,协方差和相关系数,协方差是衡量两个变量的总体误差【这个解释好像并没有什么用?】,方差是协方差的特殊情况,衡量数据集离散程度,在协方差上有这么一句话:协方差作为描述X和Y相关程度的量,在同一物理量纲之下有一定的作用,但同样的两个量采用不同的量纲使它们的协方差在数值上表现出很大的差异,为此引入相关系数,我觉得这句话对于协方差的意义更容易理解。

对DataFrame中的指定列(Series)取前k个最大的行,对应的列的索引

1
2
3
k = 10 #number of variables for heatmap
#corrmat.nlargest表示的是在SalePrice列中选出最大的前k个对应的行,也就是说这样做之后的效果是生成了10*38的DataFrame
cols = corrmat.nlargest(k, 'SalePrice')['SalePrice'].index

使用np求协方差矩阵,并用热力图显示

1
2
3
4
5
6
7
8
9
10
cm = np.corrcoef(df_train[cols].values.T)#此处之所以要转置,是因为np针对的是行,即行代表的是变量,而列代表的是值
sns.set(font_scale=1.25)#设置字体大小
"""
详细参数参考:http://seaborn.pydata.org/generated/seaborn.heatmap.html?highlight=heatmap#seaborn.heatmap
annot: 如果为true,写一个数据到每个单元上
fmt: 当增加注释时,使用的格式
annot_kws: 当annot为true时,对应值的设置
"""
hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt='.2f', annot_kws={'size': 10}, yticklabels=cols.values, xticklabels=cols.values)
plt.show()

kaggle0104
使用sns中的pairplot绘制两两变量之间的散点图,至于在下面去掉了三对强相关中其中一个变量,强相关就是多重线性问题,在广义线性回归模型中讲到过多重线性问题的解决方案,其中之一是使用惩罚项(举到岭回归,为什么岭回归可以解决多重共线性问题,其中有公式推导,这里直观的理解就是如果有个变量是多余的,也就是说没有该变量也会达到同样的效果,而惩罚项的目的就是让权重变少或者变小,这样在惩罚项的作用下必然会将该多余变量的权重降为0),也提到删除其中之一的变量,我觉得最好不要删除而用惩罚项。

1
2
3
4
#去掉3个强相关的
cols = ['SalePrice', 'OverallQual', 'GrLivArea', 'GarageCars', 'TotalBsmtSF', 'FullBath', 'YearBuilt']
sns.pairplot(df_train[cols], size = 2.5)#画对应的列之间的散点关系图,不过没搞清楚,如果是自身之间的散点图为什么是直方图
plt.show();

kaggle0105
下面这个是计算DataFrame中各列空值的占比

1
2
3
4
5
6
7
8
9
10
total = df_train.isnull().sum().sort_values(ascending=False)#计算每列中为空的个数,并排序(降序)
"""
df_train.isnull生成的是对应df_train的DataFrame,但是值为false或true;.count()计算DataFrame中每列不为空得个数
df_train.count()会忽略掉为空的元素,这样的结果就是不为空得个数
df_train.isnull().count()将空得元素转为true,这样结果就是实际行的行数
所以说下面的df_train.isnull().count()实际上可以直接用1460来替换
"""
percent = (df_train.isnull().sum()/df_train.isnull().count()).sort_values(ascending=False)#
missing_data = pd.concat([total, percent], axis=1, keys=['Total', 'Percent'])#连接,column为Total和Percent
missing_data.head(20)

删除缺失数据个数大于1的列

1
2
df_train = df_train.drop((missing_data[missing_data['Total'] > 1]).index,1)#前面的df对应值大于1的索引,1表示删除df_train的列
df_train = df_train.drop(df_train.loc[df_train['Electrical'].isnull()].index)#删除索引对应的行,axis默认为0

后面讲到数据的正态化我在其它文章中讲到过,这里就不重复了
在正态化中关于0值的处理直接不考虑感觉有点问题,而后续的文章中正态化处理使用的是log1p函数,即log(1+x),这样就没有这个问题了
关于将某列数据正态化后和价格列画散点图和之前没有正态化对比,变得不再是锥形了,这里不知道说明了什么?好像没写完一样…

数据处理及线性回归在训练集上的比较

这篇文章可以作为广义线性回归的实际篇,而且对数据处理上做的特别好
对空值的处理上做的很详细,分析每一列在添加空值,对数据集中有些分类特征列用数字表示,这样显然是不对的,因为这些数字大小并不能说明问题,这样的列是需要从数值列转换为字符串列(分类列)的

1
2
3
4
5
6
7
8
9
10
11
12
#这个做的挺好的,就是有些本来是分类特征的,却写成了数值特征,那么这是需要转换回来的
#其中key作为要修改的值,value作为替换值;(用value替换成key)
"""
转换df中的值
"""
train = train.replace({"MSSubClass" : {20 : "SC20", 30 : "SC30", 40 : "SC40", 45 : "SC45",
50 : "SC50", 60 : "SC60", 70 : "SC70", 75 : "SC75",
80 : "SC80", 85 : "SC85", 90 : "SC90", 120 : "SC120",
150 : "SC150", 160 : "SC160", 180 : "SC180", 190 : "SC190"},
"MoSold" : {1 : "Jan", 2 : "Feb", 3 : "Mar", 4 : "Apr", 5 : "May", 6 : "Jun",
7 : "Jul", 8 : "Aug", 9 : "Sep", 10 : "Oct", 11 : "Nov", 12 : "Dec"}
})

而有限类别类是可以用数值来表示的(表现在该特征的变化值),而不用最后使用独热编码(一个特征列变成多个特征了),如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
train = train.replace({"Alley" : {"Grvl" : 1, "Pave" : 2},
"BsmtCond" : {"No" : 0, "Po" : 1, "Fa" : 2, "TA" : 3, "Gd" : 4, "Ex" : 5},
"BsmtExposure" : {"No" : 0, "Mn" : 1, "Av": 2, "Gd" : 3},
"BsmtFinType1" : {"No" : 0, "Unf" : 1, "LwQ": 2, "Rec" : 3, "BLQ" : 4,
"ALQ" : 5, "GLQ" : 6},
"BsmtFinType2" : {"No" : 0, "Unf" : 1, "LwQ": 2, "Rec" : 3, "BLQ" : 4,
"ALQ" : 5, "GLQ" : 6},
"BsmtQual" : {"No" : 0, "Po" : 1, "Fa" : 2, "TA": 3, "Gd" : 4, "Ex" : 5},
"ExterCond" : {"Po" : 1, "Fa" : 2, "TA": 3, "Gd": 4, "Ex" : 5},
"ExterQual" : {"Po" : 1, "Fa" : 2, "TA": 3, "Gd": 4, "Ex" : 5},
"FireplaceQu" : {"No" : 0, "Po" : 1, "Fa" : 2, "TA" : 3, "Gd" : 4, "Ex" : 5},
"Functional" : {"Sal" : 1, "Sev" : 2, "Maj2" : 3, "Maj1" : 4, "Mod": 5,
"Min2" : 6, "Min1" : 7, "Typ" : 8},
"GarageCond" : {"No" : 0, "Po" : 1, "Fa" : 2, "TA" : 3, "Gd" : 4, "Ex" : 5},
"GarageQual" : {"No" : 0, "Po" : 1, "Fa" : 2, "TA" : 3, "Gd" : 4, "Ex" : 5},
"HeatingQC" : {"Po" : 1, "Fa" : 2, "TA" : 3, "Gd" : 4, "Ex" : 5},
"KitchenQual" : {"Po" : 1, "Fa" : 2, "TA" : 3, "Gd" : 4, "Ex" : 5},
"LandSlope" : {"Sev" : 1, "Mod" : 2, "Gtl" : 3},
"LotShape" : {"IR3" : 1, "IR2" : 2, "IR1" : 3, "Reg" : 4},
"PavedDrive" : {"N" : 0, "P" : 1, "Y" : 2},
"PoolQC" : {"No" : 0, "Fa" : 1, "TA" : 2, "Gd" : 3, "Ex" : 4},
"Street" : {"Grvl" : 1, "Pave" : 2},
"Utilities" : {"ELO" : 1, "NoSeWa" : 2, "NoSewr" : 3, "AllPub" : 4}}
)

补一下,上面是直接对df操作,为了更清楚,下面用Series:

1
all_data['Pclass']=all_data['Pclass'].replace({1:'upper',2:'middle',3:'lower'})

用下面三种方法创造新的特征(注意是增加列,而不是替换)
Simplifications of existing features
Combinations of existing features
Polynomials on the top 10 existing features
简化特征,将差别不大的归为一类

1
2
3
4
5
# 1* Simplifications of existing features
train["SimplOverallQual"] = train.OverallQual.replace({1 : 1, 2 : 1, 3 : 1, # bad
4 : 2, 5 : 2, 6 : 2, # average
7 : 3, 8 : 3, 9 : 3, 10 : 3 # good
})

合并特征

1
2
3
4
5
# 2* Combinations of existing features
# Overall quality of the house
train["OverallGrade"] = train["OverallQual"] * train["OverallCond"]
# Overall quality of the garage
train["GarageGrade"] = train["GarageQual"] * train["GarageCond"]

在和价格相关性最大的十个特征上做多项式,注意对特征多项式的方法:可以看到一个特征被开根号,平方,三次方,最后都变为一列特征

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 3* Polynomials on the top 10 existing features
train["OverallQual-s2"] = train["OverallQual"] ** 2
train["OverallQual-s3"] = train["OverallQual"] ** 3
train["OverallQual-Sq"] = np.sqrt(train["OverallQual"])
train["AllSF-2"] = train["AllSF"] ** 2
train["AllSF-3"] = train["AllSF"] ** 3
train["AllSF-Sq"] = np.sqrt(train["AllSF"])
train["AllFlrsSF-2"] = train["AllFlrsSF"] ** 2
train["AllFlrsSF-3"] = train["AllFlrsSF"] ** 3
train["AllFlrsSF-Sq"] = np.sqrt(train["AllFlrsSF"])
train["GrLivArea-2"] = train["GrLivArea"] ** 2
train["GrLivArea-3"] = train["GrLivArea"] ** 3
train["GrLivArea-Sq"] = np.sqrt(train["GrLivArea"])
train["SimplOverallQual-s2"] = train["SimplOverallQual"] ** 2
train["SimplOverallQual-s3"] = train["SimplOverallQual"] ** 3

偏度判断并正态化

1
2
3
4
5
skewness = train_num.apply(lambda x: skew(x))
skewness = skewness[abs(skewness) > 0.5]
print(str(skewness.shape[0]) + " skewed numerical features to log transform")
skewed_features = skewness.index
train_num[skewed_features] = np.log1p(train_num[skewed_features])

最后独热编码

1
train_cat = pd.get_dummies(train_cat)

交叉验证,用方均根差来计算

1
2
3
4
5
6
7
8
scorer = make_scorer(mean_squared_error, greater_is_better = False)#定义一个得分机制,均方误差,越小越好
def rmse_cv_train(model):
rmse= np.sqrt(-cross_val_score(model, X_train, y_train, scoring = scorer, cv = 10))
return(rmse)
def rmse_cv_test(model):
rmse= np.sqrt(-cross_val_score(model, X_test, y_test, scoring = scorer, cv = 10))
return(rmse)

线性回归模型

1
2
3
4
lr = LinearRegression()
lr.fit(X_train, y_train)
print("RMSE on Training set :", rmse_cv_train(lr).mean())
print("RMSE on Test set :", rmse_cv_test(lr).mean())

RMSE on Training set : 0.389629944224
RMSE on Test set : 0.4179630667
然后用岭回归模型,这里需要指定惩罚项前的系数,具体的方法是使用下面的函数

1
2
3
ridge = RidgeCV(alphas = [0.01, 0.03, 0.06, 0.1, 0.3, 0.6, 1, 3, 6, 10, 30, 60])
ridge.fit(X_train, y_train)
alpha = ridge.alpha_

该函数会在给定的列表中选取最佳的惩罚项系数,这里获得最佳惩罚项系数后继续缩小范围,如下:

1
2
3
4
5
6
ridge = RidgeCV(alphas = [alpha * .6, alpha * .65, alpha * .7, alpha * .75, alpha * .8, alpha * .85,
alpha * .9, alpha * .95, alpha, alpha * 1.05, alpha * 1.1, alpha * 1.15,
alpha * 1.25, alpha * 1.3, alpha * 1.35, alpha * 1.4],
cv = 10)#kfold=10
ridge.fit(X_train, y_train)
alpha = ridge.alpha_

岭回归模型的结果如下:

1
2
print("Ridge RMSE on Training set :", rmse_cv_train(ridge).mean())
print("Ridge RMSE on Test set :", rmse_cv_test(ridge).mean())

Ridge RMSE on Training set : 0.115405723285
Ridge RMSE on Test set : 0.116427213778
可以看到相比于线性回归模型,方均误差已经小了很多,对比训练集结果和测试机结果可以看到模型的稳定性也得到了提高(即降低了过拟合)
用Lasso回归模型,该模型使用的是L1正则化,由于L1正则化没有使用平方,而导致最后拟合的结果中很多特征权重都会降为0,这对很多无关特征的高维数据集是有效的。同上原理,得出来得结果如下:
Lasso RMSE on Training set : 0.114111508375
Lasso RMSE on Test set : 0.115832132218
这个例子说明了对于高维无关数据集上的拟合Lasso回归模型更优于岭回归模型
最后提到了elasticNet模型,这里就不说明了

xgboost+线性回归

前面的数据分析及处理部分上面都包括了,这里主要讲一下xgboost的使用,如下使用xgboost做交叉验证,图中显示的是训练集误差和测试集误差的关系,显然越往后模型过拟合越明显:

1
2
3
4
5
6
7
import xgboost as xgb
dtrain = xgb.DMatrix(X_train, label = y)
dtest = xgb.DMatrix(X_test)
params = {"max_depth":2, "eta":0.1}
#迭代次数为500,如果在100的时候误差不能继续下降就停止迭代
model = xgb.cv(params, dtrain, num_boost_round=500, early_stopping_rounds=100)#返回的是DataFrame
print(model.head())#包含训练集误差的均值和标准差,测试数据误差的均值和标准差【迭代500次后的变换】

test-rmse-mean test-rmse-std train-rmse-mean train-rmse-std
0 10.380138 0.007834 10.380139 0.003413
1 9.344810 0.008219 9.344813 0.003153
2 8.413089 0.008584 8.413092 0.002927
3 7.574992 0.008599 7.574618 0.002707
4 6.820292 0.008262 6.819935 0.002493

1
2
model.loc[:,["test-rmse-mean", "train-rmse-mean"]].plot()
plt.ylim(0,1)

kaggle0106
最后做了一下简单的模型相加,效果还不错,能够达到0.120

1
2
3
preds = 0.8*lasso_preds + 0.2*xgb_preds
solution = pd.DataFrame({"id":test.Id, "SalePrice":preds})
solution.to_csv("ridge_sol.csv", index = False)

stacking的使用

这里需要继承scikit-learn的基类,BaseEstimator为估计的基础类,RegressorMixin是回归的基类,TransformerMixin是转换类的基类
其实下面这个自定义类就是将每个模型取个平均

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
class AveragingModels(BaseEstimator, RegressorMixin, TransformerMixin):
def __init__(self, models):
self.models = models
# we define clones of the original models to fit the data in
def fit(self, X, y):
self.models_ = [clone(x) for x in self.models]
# Train cloned base models
for model in self.models_:
model.fit(X, y)
return self
#Now we do the predictions for cloned models and average them
def predict(self, X):
predictions = np.column_stack([
model.predict(X) for model in self.models_
])
return np.mean(predictions, axis=1)

使用上面的模型训练数据:

1
2
3
averaged_models = AveragingModels(models = (ENet, GBoost, KRR, lasso))
score = rmsle_cv(averaged_models)
print(" Averaged base models score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

相对于单模型提高了进一个百分点,但是上面这个模型只是简单的对每个模型进行平均就能达到这样的效果
继续使用正宗的stacking方法,如下图所示,其原理就是用交叉验证的方法来预测训练集,比如5折,进行5次【这里代码中给出的是5个相同的模型,而图中显示的是不同的模型】,每次用其他的4部分训练,取预测剩下的一部分,这样5次后就会生成一个完整的对训练集的预测结果,分别用多个模型进行上面的步骤,就会得到多组对训练集预测的结果集,最后用基模型训练这个结果集,结果集作为自变量(x),训练数据的因变量(y)作为该模型的因变量。预测的时候,先使用5次中的模型来预测测试集,并取均值;然后再用多个模型进行上述过程,得到多组预测集,再使用基模型对预测集进行预测,即为最终的结果,如下图:
kaggle0107
代码如下(此处对于5折并没有使用多模型,我觉得这样训练的模型更加稳定):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):
def __init__(self, base_models, meta_model, n_folds=5):
self.base_models = base_models
self.meta_model = meta_model
self.n_folds = n_folds
# We again fit the data on clones of the original models
def fit(self, X, y):
self.base_models_ = [list() for x in self.base_models]
self.meta_model_ = clone(self.meta_model)
kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=156)
# Train cloned base models then create out-of-fold predictions
# that are needed to train the cloned meta-model
out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))
for i, model in enumerate(self.base_models):
for train_index, holdout_index in kfold.split(X, y):
instance = clone(model)
self.base_models_[i].append(instance)
instance.fit(X[train_index], y[train_index])
y_pred = instance.predict(X[holdout_index])
out_of_fold_predictions[holdout_index, i] = y_pred
# Now train the cloned meta-model using the out-of-fold predictions as new feature
self.meta_model_.fit(out_of_fold_predictions, y)
return self
#Do the predictions of all base models on the test data and use the averaged predictions as
#meta-features for the final prediction which is done by the meta-model
def predict(self, X):
meta_features = np.column_stack([
np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)
for base_models in self.base_models_ ])
return self.meta_model_.predict(meta_features)

最后用了stacking后,将stacking的结果和Lightgbm和XGBoost的结果做了简单的加权:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
stacked_averaged_models = StackingAveragedModels(base_models = (ENet, GBoost, KRR),meta_model = lasso)
stacked_averaged_models.fit(train.values, y_train)
stacked_train_pred = stacked_averaged_models.predict(train.values)
stacked_pred = np.expm1(stacked_averaged_models.predict(test.values))
model_xgb.fit(train, y_train)
xgb_train_pred = model_xgb.predict(train)
xgb_pred = np.expm1(model_xgb.predict(test))
model_lgb.fit(train, y_train)
lgb_train_pred = model_lgb.predict(train)
lgb_pred = np.expm1(model_lgb.predict(test.values))
ensemble = stacked_pred*0.70 + xgb_pred*0.15 + lgb_pred*0.15
sub = pd.DataFrame()
sub['Id'] = test_ID
sub['SalePrice'] = ensemble
sub.to_csv('sub.csv',index=False)

确实觉得前面那篇数据处理部分做的太好了,只不过我改不过来了(太麻烦了),就只是简单的将第三部分的结果和第四部分的结果加权了一下,得出的结果由第四部分的0.1511变为0.1510

到此就先结束了。

如果觉得有帮助,给我打赏吧!