项目背景

目的：基于所提供的数据集，对房价进行预测

#首先，导入一些所需要的工具包
import numpy as np 
import pandas as pd 
%matplotlib inline
import matplotlib.pyplot as plt  
import seaborn as sns
sns.set_style('darkgrid')
import warnings
def ignore_warn(*args, **kwargs): 
    pass
warnings.warn = ignore_warn #忽略主要来自sklearn & seaborn 的警告
from scipy import stats
from scipy.stats import norm, skew # 正态分布与偏度识别
#读取数据集
train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
train.head()   #默认读取前5行

#看看测试数据，测试数据的变量和训练数据是一致的，不过会少了一列Y值
test.head()

一、数据可视化探索

1、缺失值

# id这一列对于模型预测是没有用的，我们先将它在训练集和测试集中去掉
train_ID = train['Id']
test_ID = test['Id']
train.drop("Id", axis = 1, inplace = True)     #axis=1表示列    inplace=1表示不创建新的对象，直接对象进行修改
test.drop("Id", axis = 1, inplace = True)
print("\nThe train data size after dropping Id feature is : {} ".format(train.shape)) 
print("The test data size after dropping Id feature is : {} ".format(test.shape))

#合并train test数据一起处理 
ntrain = train.shape[0]     #pd.shape[0]读取矩阵的行数
ntest = test.shape[0]
y_train = train.SalePrice.values
all_data = pd.concat((train, test)).reset_index(drop=True)   #reset_index表示重置索引
all_data.drop(['SalePrice'], axis=1, inplace=True)
print("all_data size is : {}".format(all_data.shape))             #shape读取矩阵的形状

all_data_na = (all_data.isnull().sum() / len(all_data)) * 100    #按列查看缺失率
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]   #降序查看前30个
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head(20)

#将缺失度用图表的方式展示
f, ax = plt.subplots(figsize=(15, 12))
plt.xticks(rotation='90')
sns.barplot(x=all_data_na.index, y=all_data_na)   #条形图
plt.xlabel('Features', fontsize=15)
plt.ylabel('Percent of missing values', fontsize=15)
plt.title('Percent missing data by feature', fontsize=15)

2、相关性探索

热力图是一个快速查看变量之间的相关性，帮助我们理解数据的好办法，根据下图图例，颜色越浅的地方说明两个变量的相关性越强。我们观察最后一列saleprice和其他变量的关系发现GrLivArea（地上居住面积）、OverallQual（整体质量）、GarargeCars（车库能装几辆车）这些变量的颜色比较浅，说明它们对房价的预测能力可能比较强。

#计算所有特征值，每两个之间的相关系数
corrmat = train.corr()   #列之间的相关系数
f, ax = plt.subplots(figsize=(12, 9))
sns.heatmap(corrmat, vmax=.8, square=True);   #square设置热力图矩阵小块形状，默认为False 

#选出10个与房价相关性最强的变量查看相关性系数

k = 10 #选择变量的个数
cols = corrmat.nlargest(k, 'SalePrice')['SalePrice'].index   #pandas.nlargest()输出最大相关排序，排完后选取saleprice列的索引
cm = np.corrcoef(train[cols].values.T)   #设置数据集
sns.set(font_scale=1.25)
#annot_kws设置热力图举证上数字大小字体颜色
hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt='.2f', annot_kws={'size': 10}, yticklabels=cols.values, xticklabels=cols.values)
plt.show()

通过上面的热力图，可以有如下的结论：

1. 'OverallQual', 'GrLivArea' 和 'TotalBsmtSF'与房价的相关性很强，可以后面再深入探索
2. 'GarageCars'（车库能放多少量车） 和 'GarageArea' （车库面积）和房价同样有比较强的相关性，但注意到这两个变量本身的相关性也很强，因为车库面积和车库能放多少车本身就是有强相关性的，所以这就会涉及到我们之前提到过的多重共线性的问题，对于这两个变量，我们可以去掉一个，只留一个就可以了
3. 'TotalBsmtSF' 和 '1stFloor' 似乎都是代表地下室面积，这里不确定这两个变量含义的区别，不过它们也有很强的相关性

就发现的和房价相关性较强的变量，再深入看一下它们和房价的关系

#GrLivArea代表含义是居住面积，发现和房价有比较明显的正相关关系，相关系数为0.71
#我们可以在右下角看到两个明显的outlier，价格很低，但面积巨大，我们可以考虑删除掉这两个异常值。
fig, ax = plt.subplots()
ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])  #scatter散点图
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

#TotalBsmtSF含义为地下室面积，发现地下室面积与房价似乎有更强的潜在线性关系，相关系数为0.61，同时在右侧似乎也有一个异常值存在
var = 'TotalBsmtSF'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
data.plot.scatter(x=var, y='SalePrice', ylim=(0,800000));

#GarageArea含义为车库面积，可以看到车库面积和房价也存在一定的正相关关系，相关系数为0.62
var = 'GarageArea'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
data.plot.scatter(x=var, y='SalePrice', ylim=(0,800000));

对于连续型变量，可以用散点图的方式查看它们和房价之间的关系，下面再来看一下离散型变量，对于离散型变量，可以用箱型图去查看特征不同取值中房价的分布情况

# 首先看一下房屋质量和房价之间的关系，相关系数为0.79，可以看到随着房屋整体质量的变好，房屋的整体价格也在逐渐提高，是一个比较强力的预测变量
var = 'OverallQual'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
f, ax = plt.subplots(figsize=(8, 6))
fig = sns.boxplot(x=var, y="SalePrice", data=data)  #箱线图
fig.axis(ymin=0, ymax=800000);

#再看一下房屋的建造时间和价格的关系，可以看到房屋的建造时间虽然没有和房价有明显的线性关系，但最近建造的新房屋，整体的价格会比较高，这个在后面的特征工程中会用到，注意

ar = 'YearBuilt'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
f, ax = plt.subplots(figsize=(16, 8))
fig = sns.boxplot(x=var, y="SalePrice", data=data)
fig.axis(ymin=0, ymax=800000);
plt.xticks(rotation=90);

#上面是一些比较容易想到的变量，我们可以再看一下壁炉数量和房价的关系，壁炉越多，房价越高
var = 'Fireplaces'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
f, ax = plt.subplots(figsize=(16, 8))
fig = sns.boxplot(x=var, y="SalePrice", data=data)
fig.axis(ymin=0, ymax=800000);
plt.xticks(rotation=90);

二、数据处理

通过可视化探索，我们对数据可以建立起基本的了解，下一步我们对数据的质量做一些处理。首先是异常值，我们已经发现在房屋面积和地下室面积这两个特征里面可能存在异常值，可以把它们删掉

#删除掉异常值
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index)

##对异常值进行删除重新画图，可以看到异常值没有了，在这里我们只看了一个变量，有兴趣的同学也可以看一下其他的变量异常值的分布情况。
fig, ax = plt.subplots()
ax.scatter(train['GrLivArea'], train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

#对于地下室面积，做同样处理
train = train.drop(train[(train['TotalBsmtSF']>5000) & (train['SalePrice']<200000)].index)

fig, ax = plt.subplots()
ax.scatter(train['TotalBsmtSF'], train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('TotalBsmtSF', fontsize=13)
plt.show()

下面我们来看一下缺失值的处理

#再来看一下异常值的分布情况
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]   #降序查看前30个
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head(50)

#PoolQC、MiscFeature、Alley的缺失值都在90%以上，可以考虑直接删掉这些特征
all_data= all_data.drop('PoolQC',axis=1)
all_data = all_data.drop('MiscFeature',axis=1)
all_data= all_data.drop('Alley',axis=1)
#删掉这3列后，变量变成了76列
all_data.shape

#Fence栅栏，FireplaceQu壁炉，如果数据缺失的话可能是代表房屋没有栅栏或壁炉， 我们用None填补缺失值，代表没有
all_data["Fence"] = all_data["Fence"].fillna("None")
all_data["FireplaceQu"] = all_data["FireplaceQu"].fillna("None")

#LotFrontage代表房屋前街道的长度，从现实中考虑，房屋前街道的长度可能是会和房屋所在一个街区的房屋相同，所以我们可以考虑同一个街区房屋
#LotFrontage的均值来填补缺失值
#groupby 输出所有数据中不同街区
#transform对主句进行某种同一处理（归一化、标准化等）
all_data["LotFrontage"] = all_data.groupby("Neighborhood")["LotFrontage"].transform(
    lambda x: x.fillna(x.median()))   #输出x均值填补后的x

#Garage相关的车库变量，注意到这些变量的缺失比例是完全相同的，缺失这些变量的房屋可能是没有车库，用None填充
for col in ('GarageType', 'GarageFinish', 'GarageQual', 'GarageCond'):
    all_data[col] = all_data[col].fillna('None')

#同样猜测缺失值缺失的原因可能是因为房屋没有车库，对于连续型变量用0填充
for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'):
    all_data[col] = all_data[col].fillna(0)

#地下室相关连续变量，缺失同样认为房屋可能是没有地下室，用0填充
for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath'):
    all_data[col] = all_data[col].fillna(0)

#地下室相关离散变量，同理用None填充
for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
    all_data[col] = all_data[col].fillna('None')

#Mas为砖石结构相关变量，缺失值我们同样认为是没有砖石结构，用0和none填补缺失值
all_data["MasVnrType"] = all_data["MasVnrType"].fillna("None")
all_data["MasVnrArea"] = all_data["MasVnrArea"].fillna(0)

#MSZoning代表房屋所处的用地类型，给大家先看一下不同取值的个数
#可以看到RL的取值最多
all_data.groupby('MSZoning')['MasVnrType'].count().reset_index()

业务上讲，房屋应该都有所在用地类型，且应该是上面几个值中的一个，所以有缺失值可能不是因为房屋没有用地类型导致，对于这种情况我们就最好不要用None填充，可以考虑用众数

all_data['MSZoning'] = all_data['MSZoning'].fillna(all_data['MSZoning'].mode()[0])   #mode众数

all_data = all_data.drop(['Utilities'], axis=1)
all_data["Functional"] = all_data["Functional"].fillna("Typ")
all_data['Electrical'] = all_data['Electrical'].fillna(all_data['Electrical'].mode()[0])
all_data['KitchenQual'] = all_data['KitchenQual'].fillna(all_data['KitchenQual'].mode()[0])
all_data['Exterior1st'] = all_data['Exterior1st'].fillna(all_data['Exterior1st'].mode()[0])
all_data['Exterior2nd'] = all_data['Exterior2nd'].fillna(all_data['Exterior2nd'].mode()[0])
all_data['SaleType'] = all_data['SaleType'].fillna(all_data['SaleType'].mode()[0])
all_data['MSSubClass'] = all_data['MSSubClass'].fillna("None")

#再来看看还有没有缺失值，发现没有了，到这我们就完成了对缺失值的处理
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head()

#检测因变量是否符合正态性分布
sns.distplot(train['SalePrice'] , fit=norm);  #distplot直方图

(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],  loc='best')#添加图例
plt.ylabel('Frequency')
plt.title('SalePrice distribution')

fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)   #默认检测是正态分布
plt.show()

从房价的直方图分布和qq-plot可以看出，房价的分布是右偏的，需要对其做一些转换让它符合正太分布

采用log对数变换对房价进行处理，通过下面转换后房价的分布可以看出，房价转换后符合正态分布

train["SalePrice"] = np.log1p(train["SalePrice"])    #数据平滑处理，使其更加服从高斯分布

sns.distplot(train['SalePrice'] , fit=norm);    #distplot直方图

(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)], loc='best')#添加图例
plt.ylabel('Frequency')
plt.title('SalePrice distribution')

fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)   #默认检测是正态分布
plt.show()

三、特征工程

对数据做完基本的数据处理后，下面就进入我们最重要最重要的特征工程！

#1、特征值处理  2、特征提取  3、特征选择
#下面我们按照上面的3个步骤，开始进行特征工程，首先我们尝试创造出一些新特征,用下面两个例子给大家做个示范
#加地下室面积、1楼面积、2楼面积加起来可以得到房屋总面积特征
all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF']
#之前可视化时候注意到，建造时间比较近的房子房价比较高，所以新创造一个01特征，如果房屋建造时间在1990年后，则为1，否则是0
all_data['YearBuilt_cut'] = all_data['YearBuilt'].apply(lambda x:1 if x>1990 else 0)

# 我们看一下这两个新加的特征与price的关系
tep=all_data[:ntrain]
tep['SalePrice']=y_train

#房屋总面积特征，发现与房价有比较强的线性关系，但也有2个可能的异常值，可以考虑再去掉这两个异常值
fig, ax = plt.subplots()
ax.scatter(tep['TotalSF'], tep['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('TotalSF', fontsize=13)
plt.show()

#建筑年限，可以看到1990年前建造的房子和1990年后建造的房子房价再分布上有较大的差异
var='YearBuilt_cut'
data=pd.concat([tep['SalePrice'], tep[var]], axis=1)
f,ax=plt.subplots(figsize=(16,8))
fig=sns.boxplot(x=var, y='SalePrice', data=data)
fig.axis(ymin=0, ymax=800000)
plt.xticks(rotation=90)

#对有序离散变量，用label encoder进行编码
from sklearn.preprocessing import LabelEncoder
cols=("FireplaceQu", "BsmtQual", "BsmtCond", "GarageQual", "GarageCond", 'ExterQual', "ExterCond", "HeatingQC", "KitchenQual", 'BsmtFinType1',
     'BsmtFinType2', 'Functional', 'BsmtExposure', 'GarageFinish', 'LandSlope','LotShape', 'PavedDrive', 'Street', 'CentralAir', 'MSSubClass', 'OverallCond', 
        'YrSold', 'MoSold')

# process columns, apply LabelEncoder to categorical features
for c in cols:
    lbl=LabelEncoder()
    lbl.fit(list(all_data[c].values))
    all_data[c]=lbl.transform(list(all_data[c].values))
    
#shape
print('Shape all_data: {}'.format(all_data.shape))

#再来看一下数据中的格式，发现之前的英文字符串已经变成了数字
all_data.head()

#数据集中还有部分非有序性离散变量，我们将他们转换成哑变量的形式（和onehot一个意思）
all_data=pd.get_dummies(all_data)
all_data.shape

all_data.shape

下面进行特征工程的最后一步，特征筛选，为避免多重共线性问题，下面我们会识别皮尔森相关性系数大于0.9的特征，并将这些特征删除。

threshold=0.9
#相关性矩阵
corr_matrix=all_data.corr().abs()
corr_matrix.head()

#只选择矩阵的上半部分
#np.triu中k表示从第几条对角线起保留数据，k=0表示主对角线
#df.where从主体df出发，true返回df本身
upper=corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))
upper.head()

#有6列特征需要删除
to_drop=[column for column in upper.columns if any(upper[column]>threshold)]

print('There are %d columns to remove.' % (len(to_drop)))

all_data=all_data.drop(columns=to_drop)
all_data.shape

#之前是把train和test数据集放在一起处理，现在再把他们分开
train=all_data[:ntrain]
test=all_data[ntrain:]

四、建模

在做模型之前，为了衡量模型的泛化误差避免出现过拟合等现象，先要建立一个本地验证集。建模是用训练数据训练模型，并且期望模型在推广到其他数据上时，也有比较好的预测效果，在比赛中对应的就是测试集。但由于测试集没有Y值，即没有真实的房价，所以我们没有办法对模型在测试集上的预测准确率做出判断。因此接下来将训练集八二分，留出20%的数据去验证模型的效果，作为验证集来测试模型的泛化效果。

#导入一些所需要的数据包
from sklearn.linear_model import Ridge, RidgeCV, ElasticNet, LassoCV, LassoLarsCV
from sklearn.model_selection import cross_val_score
#使用cross-val-score函数去作为模型准确性的评估
#np.sqrt()计算各元素的平方根     cv交叉验证折数或者可迭代次数
def rmse_cv(model):
    rmse=np.sqrt(-cross_val_score(model, train, y_train, scoring="neg_mean_squared_error", cv=5))
    return(rmse)
#导入ridge模型
model_ridge=Ridge()

在导入模型后，需要对模型进行调参，对于lasso模型来说，主要需要调节的参数是alpha，alpha决定模型的正则化程度。正则化程度越高，我们的模型越不容易过度拟合。然而，它也会失去灵活性，可能会欠拟合，我们需要在过拟合和欠拟合中寻找一个平衡，找到那个最佳参数。因此，可以列出不同的alpha取值，用上面定义的rmse_cv模型评分参数去看每个参数下模型的分数怎么样

alphas=[0.05, 0.1, 0.3, 1, 3, 5, 10, 15, 30, 50, 75]
cv_ridge=[rmse_cv(Ridge(alpha=alpha)).mean() for alpha in alphas]
cv_ridge=pd.Series(cv_ridge, index=alphas)
cv_ridge.plot(title='Validation - Just Do It')
plt.xlabel('alpha')
plt.ylabel('rmse')

可以看到曲线是一个耐克的标志，先降低后升高，在alpha值为5时，模型的误差最小。

cv_ridge

#alpha参数用我们之前验证过的5，然后用训练集对模型进行训练
clf=Ridge(alpha=5)
clf.fit(train, y_train)
#最后对测试集进行预测，我们就完成了本次项目
predict=clf.predict(test)
sub=pd.DataFrame()
sub['Id']=test_ID
sub['SalePrice']=predict
sub.head(10)