目录
1 项目背景
2 初始数据分析
目标值分析
特征与目标值相关性
变量特征相关性
3 数据预处理
目标变量正态分布化
异常值处理
缺失值处理
转换特征
保存训练集和测试集
4 模型预测
岭回归
lasso
随机森林
5 预测结果
1 项目背景
项目链接:House Prices – Advanced Regression Techniques | Kaggle
这是kaggle的一个经典Data Science项目,作为数据分析的新手,房价预测是一个很好的入门练习项目。
数据集分为训练集‘train.csv’和测试集‘test.csv’,要求根据房子的质量、面积、街区、壁炉个数等79个特征,预测相应的房价。
评价指标是回归问题中常用的均方误差(RMSE):
2 初始数据分析
首先导入所需的库:
import warningswarnings.filterwarnings("ignore") # 忽略警告信息import numpy as npimport pandas as pd from scipy.stats import norm, skew #获取统计信息from scipy import stats # 离散统计分布以及连续统计分布import seaborn as sns # 绘图包color = sns.color_palette()sns.set_style('darkgrid')import matplotlibimport matplotlib.pyplot as plt pd.set_option('display.float_format', lambda x: '{:.3f}'.format(x)) # 限制浮点输出到小数点后3位import os
读取训练集和测试集,检查数据集大小,并查看前5行数据:
# 加载数据train = pd.read_csv('kagglehouse/train.csv')test = pd.read_csv('kagglehouse/test.csv')# 检查样本和特征的数量print("训练集初始大小: {} ".format(train.shape))print("测试集初始大小: {} ".format(test.shape))
训练集初始大小: (1460, 81)
测试集初始大小: (1459, 80)
#查看前5行数据print(train.head())
print(test.head())
训练集初始大小为1460行81列,测试集初始大小为1459行80列,训练集比测试集多出最后一列价格列,这也是我们测试集的预测目标。
其中第一列Id列没有用,直接删掉:
# 删除原数据集的Id列train.drop("Id", axis=1, inplace=True)test.drop("Id", axis=1, inplace=True)
目标值分析
目标值就是房屋价格,我们先看一看它的分布情况:
#绘制目标值分布sns.distplot(train['SalePrice'])plt.show()
是明显的右偏分布,和正态分布有较大差距,我们再看一下这个分布曲线的峰度和偏度:
# 计算房价的峰度和偏度SP_skew = train['SalePrice'].skew()SP_kurt = train['SalePrice'].kurt()print('峰度:',SP_skew)print('偏度:',SP_kurt)
峰度: 1.8828757597682129
偏度: 6.536281860064529
峰度(Kurtosis)是描述总体中所有取值分布形态陡缓程度的统计量,是与正态分布作比较的。
峰度为0表示该总体数据分布与正态分布的陡缓程度相同;峰度大于0表示该总体数据分布与正态分布相比较为陡峭,为尖顶峰;峰度小于0表示该总体数据分布与正态分布相比较为平坦,为平顶峰。
偏度(Skewness)描述的是某总体取值分布的对称性,同样与正态分布相比较。
偏度为0表示其数据分布形态与正态分布的偏斜程度相同;偏度大于0表示其数据分布形态与正态分布相比为正偏或右偏,即有一条长尾巴拖在右边,数据右端有较多的极端值;偏度小于0表示其数据分布形态与正态分布相比为负偏或左偏。
根据图可知,数据目标值峰度大于0,为尖顶峰,偏度较大,整体呈右偏分布,后期需要我们做一些处理,使其向正态分布靠拢,因为回归模型在正态分布的数据集上表现更好。
查看目标值统计信息:
#查看目标值统计信息print(train['SalePrice'].describe())
count 1460.000
mean 180921.196
std 79442.503
min 34900.000
25% 129975.000
50% 163000.000
75% 214000.000
max 755000.000
Name: SalePrice, dtype: float64
最大值太大了,与75%的取值都差很多,推断最大值附近应该有异常值存在。
特征与目标值相关性
下面看一下各个特征与目标值之间的关系,这里要先分开类别特征和数字特征:
num_features = []cate_features = []for col in test.columns: if test[col].dtype == 'object': cate_features.append(col) else: num_features.append(col)print('数值型特征:', len(num_features))print('类别型特征:', len(cate_features))
数值型特征: 36
类别型特征: 43
先查看各数值型特征与目标值之间的关系,采用散点图:
#数值特征与目标值画散点图plt.figure(figsize=(20, 30))plt.subplots_adjust(left=0.1, bottom=0.1,wspace=0.6,hspace=0.6)for i, feature in enumerate(num_features): plt.subplot(6,6, i+1) sns.scatterplot(x=feature, y='SalePrice', data=train, alpha=0.5) plt.xlabel(feature) plt.ylabel('SalePrice')plt.show()
通过散点图发现,“GrLivArea”(地上生活面积)、“TotalBsmtSF”(地下室总面积)、“1stFlrSF”(一层面积)、“2ndFlrSF”(二层面积)这几个数据与目标值呈现明显的线性关系,他们都是与面积相关的,后续可以构造新特征,这也符合常理。
然后再查看一些类别型变量与目标值的关系,用箱型图表示:
上图中的“OverallQual”(房子材料与装修水平)虽然是数值型数据,但是数值实际代表类别,与目标值有明显相关性:
plt.figure(figsize=(20, 16))sns.boxplot(x='OverallQual', y='SalePrice', data=train)plt.xlabel('OverallQual')plt.ylabel('SalePrice')plt.xticks(rotation=90, fontsize=12)plt.show()
可见装修水平越高,房子售价越高。但是同时意味着购买的人越多么?显然不是,我们通过计算不同装修水平的售出价总和判断哪个装修水平最受青睐:
#柱状图OverallQual_SalePrice = train[['OverallQual','SalePrice']].groupby(['OverallQual']).sum().sort_values(by='OverallQual',ascending=True)OverallQual_SalePrice.plot(kind='bar')plt.show()
可见5、6、7级装修水平最受买房者青睐。
下面看一下重要的“Neighborhood”(街区位置)与房屋价格的关系:
可见不同街区对应不同的房屋价格,NAmes和CollgCr地段房子销售额最高。
变量特征相关性
接下来看一下变量之间的相关性:
#多变量分析corrs = train.corr()plt.figure(figsize=(20, 9))sns.heatmap(corrs)plt.show()
由于特征较多,选取与目标值相关度最高的十个变量再次分析:
corrs = train.corr()cols_10 = corrs.nlargest(10, 'SalePrice')['SalePrice'].indexcorrs_10 = train[cols_10].corr()plt.figure(figsize=(6, 6))sns.heatmap(corrs_10, annot=True)plt.show()
3 数据预处理
目标变量正态分布化
由前面分析已经看到,目标变量呈右偏分布,我们要把他调到标准正态分布,可以采用对数变换。线性回归模型要求变量为正态分布。
y = train['SalePrice']y = np.log1p(y)sns.distplot(y, fit=norm)print('处理后峰度:',y.skew())print('处理后偏度:',y.kurtosis())plt.show()
处理后峰度: 0.12134661989685333
处理后偏度: 0.809519155707878
异常值处理
#处理异常值sns.scatterplot(x='GrLivArea', y='SalePrice', data=train)plt.show()
发现右下角有明显的异常值,我们去掉这个异常值:
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<200000)].index)sns.scatterplot(x='GrLivArea', y='SalePrice', data=train)plt.show()
同时发现TotalBsmtSF中的异常值也被去掉了:
缺失值处理
很多人是将训练集和测试集拼在一起进行处理,但是这样实际上造成了数据泄露,应该先处理训练值,再用训练值数据同样处理测试值。
先看一下训练集的缺失值个数:
train_nan = train.isnull().sum()train_nan = train_nan.drop(train_nan[train_nan==0].index).sort_values(ascending=False)print(train_nan)
测试集情况也差不多。下面我们填补缺失值。
类别特征缺失值
在data_description.txt中已有说明,一部分特征值的缺失是因为这些房子根本没有该项特征,对于这种情况我们统一用“None”或者“0”来填充。
Python中的None是一个特殊常量,不是0,也不是False,不是空字符串,更多的是表示一种不存在,是真正的空。它就是一个空的对象,只是没有赋值而已。如PoolQC确实的意思是表示“没有泳池”。
以下这些类别的缺失值填充为None:
none_lists = ['PoolQC', 'MiscFeature', 'Alley', 'Fence', 'FireplaceQu', 'GarageType', 'GarageFinish', 'GarageQual', 'GarageCond', 'BsmtFinType1', 'BsmtFinType2', 'BsmtCond', 'BsmtExposure', 'BsmtQual', 'MasVnrType']for col in none_lists: train[col] = train[col].fillna('None') test[col] = test[col].fillna('None')
对于缺失较少的离散型特征,比如Electrical,用众数填补缺失值:
#填充众数most_lists = ['MSZoning', 'Exterior1st', 'Exterior2nd', 'SaleType', 'KitchenQual', 'Electrical']for col in most_lists: train[col] = train[col].fillna(train[col].mode()[0]) test[col] = test[col].fillna(train[col].mode()[0])
Functional居家功能性,数据描述说NA表示类型”Typ”。因此,我们用其填充。
#Functional填充Typtrain['Functional'] = train['Functional'].fillna('Typ')test['Functional'] = test['Functional'].fillna('Typ')
数值特征缺失值
填充0
zero_lists = ['GarageYrBlt', 'MasVnrArea', 'BsmtFullBath', 'BsmtHalfBath', 'BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF', 'GarageCars', 'GarageArea', 'TotalBsmtSF']for col in zero_lists: train[col] = train[col].fillna(0) test[col] = test[col].fillna(0)
LotFrontage,由于房屋到街道的距离,最有可能与其附近其他房屋到街道的距离相同或相似,因此我们可以通过该社区的LotFrontage中位数来填充缺失值:
train['LotFrontage'] = train.groupby('Neighborhood')['LotFrontage'].apply(lambda x: x.fillna(x.median()))for ind in test['LotFrontage'][test['LotFrontage'].isnull().values==True].index: x = test['Neighborhood'].iloc[ind] test['LotFrontage'].iloc[ind] = train.groupby('Neighborhood')['LotFrontage'].median()[x]
最后我们删除Utilities这个特征,这个特征除了一个“NoSeWa”和 2 个 NA ,其余值都是“AllPub”,因此该项特征的方差非常小。这个特征对预测建模没有帮助。因此,我们可以安全地删除它。
train = train.drop(['Utilities'], axis=1)test = test.drop(['Utilities'], axis=1)
最后检查是否还有缺失值:
print(train.isnull().sum().max())print(test.isnull().sum().max())
返回都为0
转换特征
转换数值特征为类别特征
部分数值特征虽然是数值,但实际上表征的是类别,因此我们将其转化为类别特征:
train['MSSubClass'] = train['MSSubClass'].apply(str) # apply()函数默认对列进行操作train['YrSold'] = train['YrSold'].astype(str)train['MoSold'] = train['MoSold'].astype(str)test['MSSubClass'] = test['MSSubClass'].apply(str) test['YrSold'] = test['YrSold'].astype(str)test['MoSold'] = test['MoSold'].astype(str)
转换类别特征为数值特征
有些类别特征实际上有高低好坏之分,这些特征的质量越高,就可能在一定程度导致房价越高,因此要把他转换成编码形式。对于各个类别中可能存在顺序关系的,用LabelEncoder编码,对于不存在顺序关系的,用get_dummies。
先进行LabelEncoder编码,独热编码最后处理:
#LabelEncoder编码for col in cate_features: train[col] = train[col].astype(str) test[col] = test[col].astype(str)cols = ['Street', 'Alley', 'LotShape', 'LandContour', 'LandSlope', 'HouseStyle', 'RoofMatl', 'Exterior1st', 'Exterior2nd', 'ExterQual', 'ExterCond', 'Foundation', 'BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2', 'HeatingQC', 'CentralAir', 'KitchenQual', 'Functional', 'FireplaceQu', 'GarageFinish', 'GarageQual', 'GarageCond', 'PavedDrive', 'PoolQC', 'Fence']for col in cols: encoder = LabelEncoder() value_train = set(train[col].unique()) value_test = set(test[col].unique()) value_list = list(value_train | value_test) encoder.fit(value_list) train[col] = encoder.transform(train[col]) test[col] = encoder.transform(test[col])
偏斜特征
对于数值型特征,我们希望它们尽量服从正态分布,也就是不希望这些特征出现正负偏态。那么我们先来计算一下各个特征的偏度:
# 计算特征的偏度numeric_data = train[num_features]skewed_feats = numeric_data.apply(lambda x: skew(x.dropna())).sort_values(ascending=False)skewness = pd.DataFrame({'Skew' :skewed_feats})print(skewness.head(10))
对高偏度的特征进行Box-Cox变换:
new_skewness = skewness[skewness.abs() > 0.5]print("有{}个高偏度特征被Box-Cox变换".format(new_skewness.shape[0]))
有29个高偏度特征被Box-Cox变换
将偏度平滑至0.15
from scipy.special import boxcox1pskewed_features = new_skewness.index lam = 0.15for feat in skewed_features: train[feat] = boxcox1p(train[feat], lam) test[feat] = boxcox1p(test[feat], lam)
构建新特征
先把训练集和测试集合并,然后构造几个与目标特征相关的新特征
#把所有数据连上all_data = pd.concat((train.drop('SalePrice', axis=1), test)).reset_index(drop=True)#构造新特征all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF'] # 房屋总面积all_data['OverallQual_TotalSF'] = all_data['OverallQual'] * all_data['TotalSF'] # 整体质量与房屋总面积交互项all_data['OverallQual_GrLivArea'] = all_data['OverallQual'] * all_data['GrLivArea'] # 整体质量与地上总房间数交互项all_data['OverallQual_TotRmsAbvGrd'] = all_data['OverallQual'] * all_data['TotRmsAbvGrd'] # 整体质量与地上生活面积交互项all_data['GarageArea_YearBuilt'] = all_data['GarageArea'] * all_data['YearBuilt'] # 车库面积与建造时间交互项all_data['IsRemod'] = 1 all_data['IsRemod'].loc[all_data['YearBuilt']==all_data['YearRemodAdd']] = 0 #是否翻新(翻新:1, 未翻新:0)all_data['BltRemodDiff'] = all_data['YearRemodAdd'] - all_data['YearBuilt'] #翻新与建造的时间差(年)
独热编码
把剩余的类别特征赋予独热编码:
dummy_features = list(set(cate_features).difference(set(cols)))all_data = pd.get_dummies(all_data, drop_first=True)
保存训练集和测试集
trainset = all_data[:1458]y = train['SalePrice']trainset['SalePrice'] = y.valuestestset = all_data[1458:]trainset.to_csv('train_data.csv', index=False)testset.to_csv('test_data.csv', index=False)print('The shape of training data:', trainset.shape)print('The shape of testing data:', testset.shape)
The shape of training data: (1458, 163)
The shape of testing data: (1459, 162)
可以看到,通过一系列处理之后,训练集的列数增加到163,测试集为162。
4 模型预测
导入所需库
import numpy as npimport pandas as pdimport timeimport seaborn as snsimport matplotlib.pyplot as pltfrom sklearn.linear_model import Lasso, Ridgefrom sklearn.ensemble import RandomForestRegressorfrom sklearn.model_selection import KFold, cross_val_scorefrom sklearn.metrics import mean_squared_errorimport warningsdef ignore_warn(*args, **kwargs): passwarnings.warn = ignore_warnpd.set_option('display.float_format', lambda x: '{:.3f}'.format(x)) # 限制浮点输出到小数点后3位
岭回归
拟定不同的50个参数,通过交叉验证选取适合的参数。交叉验证,存储不同的alpha下,均方误差,通过绘图不同参数下的误差曲线查看最好的参数。
y_train = np.log1p(train.pop('SalePrice'))#岭回归alphas =np.logspace(-3,2,50)test_scores = []for alpha in alphas: clf = Ridge(alpha) test_score = np.sqrt(-cross_val_score(clf,train.values,y_train,cv=10,scoring='neg_mean_squared_error')) test_scores.append(np.mean(test_score))plt.plot(alphas,test_scores)plt.show()
进一步找到最佳位置:
alphas =np.logspace(0,1,50)
alpha=2.0附近,平均test_score最低,为0.1127
lasso
lasso_alpha = np.logspace(-4,1,50)test_scores = []for alpha in lasso_alpha: lasso =Lasso(alpha) test_score = np.sqrt(-cross_val_score(lasso, train.values, y_train, scoring='neg_mean_squared_error', cv=10)) test_scores.append(np.mean(test_score)) print(alpha,np.mean(test_score))plt.plot(lasso_alpha,test_scores)plt.show()
lasso_alpha = np.logspace(-5,-2,50)
alpha=0.0001附近,平均test_score最低,为0.1116
随机森林
主要参数:
n_estimators:表示森林里树的个数;
max_features:随机选择特征集合的子集合,并用来分割节点。
将树的个数作为变量,以子集0.3倍的特征,采用交叉验证的方式,来探索多少树个数能得到好的效果。
from sklearn.ensemble import RandomForestRegressorNs = [50,100,150,200,250,300,350,400]test_scores = []for N in Ns: clf = RandomForestRegressor(n_estimators=N, max_features=0.3) test_score = np.sqrt(-cross_val_score(clf, train.values, y_train, cv=10, scoring='neg_mean_squared_error')) test_scores.append(np.mean(test_score)) print(f"N:{N},平均test_score:{sum(test_score)/10}")plt.plot(Ns,test_scores)plt.show()
N=350时,平均test_score最低,为0.136
5 预测结果
最终选择错误率最低的lasso模型预测测试集:
lasso = Lasso(alpha=0.0001)lasso.fit(train.values, y_train)lasso_predict = lasso.predict(test.values)y_lasso = np.expm1(lasso_predict)
保存结果
test0 = pd.read_csv('kagglehouse/test.csv')test_id = test0['Id']submission=pd.DataFrame()submission['Id']=test_idsubmission['SalePrice'] = y_lassosubmission.to_csv('kagglehouse/submission.csv', index=False)
我在kaggle上上传了最终结果,排名891,大概在21%,作为第一次kaggle比赛,接下来还需要后续的优化过程~
参考资料
Kaggle经典项目——房价预测_Harold_Ran的博客-CSDN博客Kaggle房价预测详解_夜的乄第七章的博客-CSDN博客_kaggle房价预测
kaggle 房价预测_xwyzsn的博客-CSDN博客