加利福尼亚州房价分析与预测

本文通过分析加利福尼亚住房数据集，建立了一个模型来估计加利福尼亚地区的房价。

一、准备工作

1.导入必要的包

# 导入必要的包
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

2. 数据说明：

该数据集是加利福尼亚住房数据集的修改版本，可从LuísTorgo的页面（波尔图大学）获得。LuísTorgo从StatLib存储库（现已关闭）获取它。也可以从StatLib镜像下载数据集。

该数据集出现在1997年由Pace，R.Kelley和Ronald Barry撰写的名为Sparse Spatial Autoregressions的论文中，该论文发表于“ 统计与概率快报 ”期刊。他们使用1990年加州人口普查数据建立了它。每个人口普查区块组包含一行。区块组是美国人口普查局发布样本数据的最小地理单位（区块组通常拥有600至3,000人口）。

3. 下载与载入数据

你可以直接使用下面的代码下载数据，也可以访问 https://github.com/ageron/handson-ml/tree/master/datasets/housing 下载数据集

import os
import tarfile
from six.moves import urllib
DOWNLOAD_ROOT = "https://raw.githubusercontent.com/ageron/handson-ml/master/"
HOUSING_PATH = os.path.join("datasets", "housing")
HOUSING_URL = DOWNLOAD_ROOT + "datasets/housing/housing.tgz"
def fetch_data(housing_url = HOUSING_URL,housing_path = HOUSING_PATH):
    if not os.path.isdir(housing_path):
        os.makedirs(housing_path)
    tgz_path = os.path.join(housing_path,"housing_tgz")
    urllib.request.urlretrieve(housing_url,tgz_path)
    housing_tgz = tarfile.open(tgz_path)
    housing_tgz.extractall(path=housing_path)
    housing_tgz.close()

# 载入数据
data = pd.read_csv('./housing.csv')
data.sample(5)

二、数据初探 与 描述性统计分析

在这部分中，我们不会修改任何数据，只是通过各种统计手段帮助我们对这组数据有一个直观的了解，为之后的处理提取信息。

data.info()

从上面的结果中可以看到这组数据含有 20640 组实例，每个实例含有 10 个特征，包括：地理位置信息，房屋年限、总房间大小、总卧室大小、该地区人口、房屋使用者人数、平均人数、到海边的距离。其中 total_bedrooms 含有缺失值，其他的数据还是很完整的。

首先我们直观的感受一下地理位置，基本上全州都被覆盖了，有一些地方的房子很密集。

从上面的分析我们还可以注意到，大多数的特征都是数值型的变量，只有一个是字符串型的类别变量，我们先查看一下这个变量

data['ocean_proximity'].value_counts()

我们可以看到这个变量分为五个类别，描述的是到海边距离的远近，大多数都是靠海或者一小时车程，只有 5 个小区是处在内陆。

我们接下来在看看数值型的变量

data.describe()

看这个数据就要考验数据敏感性啦，我也看不出来多少，就是感觉各个特征之间的数量级差的还是很大，很多变量最大值与最小值差距也很大，这都可能会对我们之后的建模造成一定的影响。

# 绘制每个属性的直方图，来快速了解数据。
%matplotlib inline
import matplotlib.pyplot as plt

data.hist(bins=50,figsize=(20,15))
plt.show()

从这个统计图上来看还是可以看出来一些问题的，比如说 housing_median_age 和 median_house_value 超过 50 的明显高于之前的其他种类，这说明在数据集中，高于 50 的被划分为了 50，这可能会给我们带来一些困扰。另外地理位置的信息似乎并不近似满足正态分布（因为有两个峰），这也可能会影响我们最后模型的精确度。

三、创建测试集与训练集

1. 方法1 - 基于随机数的选取

第一种创建测试集与训练集的方法是通过随机指定下标创建，下面是一个使用 numpy 实现的例子，但是在实际操作中我们一般直接使用 sklearn 中的一个完美替代。

# 仅供示范
def split_train_test(data,test_ratio):
    '''
        生成训练集和测试集
    '''
    np.random.seed(42)
    #生成不重复的下标随机数
    shuffled_indices = np.random.permutation(len(data))
    test_set_size = int(len(data) * test_ratio)
    # 测试集的下标集合
    test_indices = shuffled_indices[:test_set_size]
    # 训练集的下标集合
    train_indices = shuffled_indices[test_set_size:]
    
    return data.iloc[train_indices] , data.iloc[test_indices]

train_set,test_set = split_train_test(data, 0.2)
len(train_set),len(test_set)

#(16512, 4128)

# 实际操作使用这个
from sklearn.model_selection import train_test_split
train_set , test_set = train_test_split(data, test_size=0.2, random_state=42)
len(train_set), len(test_set)

#(16512, 4128)

有的时候我们需要将标签也分开，可以使用这个函数：

X_train,y_train,X_text,y_test = train_test_split(data,target,test_size= 0.2,random_state=44)

2. 方法2 - 基于唯一标识的选取

上面的这种方法实现很容易，使用起来也很简单，但是有一个问题就是一旦我们新添加了数据，那么使用这种方法就可能将原来的测试集数据暴漏给模型训练，那么我们的测试集就没有意义了。

我们有一种想法就是选择一个唯一标识，只要这个标识不变那么选择状态就不变，最简单的方法就是根据 id % 5 == 1 选取

下面的这种方法是基于哈希算法，选取一个唯一标识进行映射，然后取映射结果的最后一个字节，这个字节大小不会超过 256。

这个方法虽然解决了上面这个问题，但是这个方法最后得到的测试集数量并不是十分准确的。

import hashlib
# 选择测试集，方式2.无惧新增数据
def test_set_check(identifier, test_radio, hash):
    return hash(np.int64(identifier)).digest()[-1] < (256 * test_radio)

def split_train_test_by_id(data, test_radio, id_column, hash=hashlib.md5):
    ids = data[id_column]
    in_test_set = ids.apply(lambda id_:test_set_check(id_,test_radio, hash))
    return data.loc[~in_test_set], data.loc[in_test_set]

# 为数据添加index列
data_width_id = data.reset_index()

train_set, test_set = split_train_test_by_id(data_width_id,0.2,"index")

len(train_set),len(test_set)

(16362, 4278)

分层随机选取

在后面的特征工程中我们会发现，平均收入（）特征对房价的影响十分巨大，这个变量大约可以解释 70% 的房价，那么我们在选取数据的时候就要保证数据集中的每个分层都要有足够的实例位于我们的训练集中。否则，最后模型的评估就会有偏差。这意味着，你不能有过多的分层，且每个分层都要足够大。后面的代码通过将收入中位数除以 1.5（以限制收入分类的数量），创建了一个收入类别属性，用ceil对值舍入（以产生离散的分类），然后将所有大于 5的分类归入到分类 5

data['income_category'] = np.ceil(data['median_income'] / 1.5) # 计算大于等于改值的最小整数
data['income_category'].where(data['income_category'] < 5 , 5.0, inplace=True)
data['income_category'].hist(bins=5,histtype='stepfilled')
plt.show()

直观上来看我们设置的参数还好，基本可以分开这五类，我们还是用数据计算一下

# 查看每个类别在总数总所占比例。
data['income_category'].value_counts() / len(data)

from sklearn.model_selection import StratifiedShuffleSplit

# 分层随机抽取测试集
split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)
result = split.split(data, data['income_category'])
for train_index,test_index in result:
    strat_train_set = data.loc[train_index]
    strat_test_set = data.loc[test_index]

len(strat_train_set) ,len(strat_test_set),

#(16512, 4128)

我们再来看一下抽取之后各个类别的比例，从结果上来看与之前的比例相近。

# 测试集中每个类别的数据比例
strat_test_set['income_category'].value_counts()/len(strat_test_set)

# 删除收入类别列
for set in (strat_test_set,strat_train_set):
    # inplace = True 时，在源数据上替换，并返回None 
    set.drop(['income_category'],axis=1,inplace=True)
strat_train_set.head()

四、特征选择与新特征的构建

第四、第五部分都是属于特征工程的内容，有一句流传了很久的话就是，特征工程的高度决定了模型的精确度，可见其重要性。但是奈何博主平时专业课没好好学，这方面所知甚少，只好使用一些极为简单和基础的办法来处理。特征选择的话我们就只是基于简单的相关性来判断。首先我们将数据可视化直观的观察一下房价大概会和什么因素有关。

housing = strat_train_set.copy()

# 将地理数据可视化
housing.plot(kind='scatter',x='longitude',y='latitude',alpha=0.4,
             s=housing['population'] / 100,label='population',c="median_house_value",
            cmap=plt.get_cmap("jet"),colorbar=True)
plt.legend()
plt.show()

上图中每一个点代表一个小区，点的面积代表人口，面积越大人口越多，颜色代表房价，红色代表房价较高，蓝色代表房价较低。从上图我们大概可以感受到，在纬度为 37 34 附近的沿海房价较高，此外人口适中的小区房价更高，沿海普遍要比内陆房价高。

下面我们计算一下相关性矩阵

mask = np.zeros_like(housing.corr(), dtype=np.bool)
#mask[np.triu_indices_from(mask)] = True

plt.subplots(figsize = (15,12))
sns.heatmap(housing.corr(), 
            annot=True,
            #mask = mask,
            cmap = 'RdBu_r',
            linewidths=0.1, 
            linecolor='white',
            vmax = .9,
            square=True)
plt.title("Correlations Among Features", y = 1.03,fontsize = 17);
plt.show()

corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)

从上面这个相关性矩阵可以看出很多问题，比如与我们所关心的房价相关性较强的变量很少，另外从上面的图中也可以看到中间有很大一块是枣红色，这意味着变量之间的相关性很强。那么我们需要对变量进行转化。

housing["rooms_per_household"] = housing['total_rooms'] / housing["households"]
housing["bedrooms_per_room"] =  housing['total_bedrooms'] / housing["total_rooms"]
housing["population_per_household"] = housing['population'] / housing["households"]

我们再看一下相关性

corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)
corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)

from sklearn.base import BaseEstimator,TransformerMixin
rooms_ix, bedrooms_ix,population_ix, househould_ix = 3, 4, 5, 6

class CombineAttributesAdder(BaseEstimator,TransformerMixin):
    def __init__(self,add_bedrooms_per_room = True):
        self.add_bedrooms_per_room = add_bedrooms_per_room
        
    def fit(self, X, Y = None):
        return self
    
    def transform(self, X, Y = None):
        rooms_per_household = X[:,rooms_ix] / X[:,househould_ix]
        population_per_household = X[:,population_ix] / X[:,househould_ix]

        if self.add_bedrooms_per_room:
            bedrooms_per_room = X[:,bedrooms_ix] / X[:,rooms_ix]
            # 将数据按列拼接
            return np.c_[X, rooms_per_household, population_per_household, bedrooms_per_room]
        else :
            return np.c_[X,rooms_per_household, population_per_household]
        

attr_adder = CombineAttributesAdder(add_bedrooms_per_room=False)
housing_extra_attribs = attr_adder.transform(housing.values)

housing_extra_attribs

array([[-121.89, 37.29, 38.0, ..., 2.094395280235988, 4.625368731563422,
        2.094395280235988],
       [-121.93, 37.05, 14.0, ..., 2.7079646017699117, 6.008849557522124,
        2.7079646017699117],
       [-117.2, 32.77, 31.0, ..., 2.0259740259740258, 4.225108225108225,
        2.0259740259740258],
       ...,
       [-116.4, 34.09, 9.0, ..., 2.742483660130719, 6.34640522875817,
        2.742483660130719],
       [-118.01, 33.82, 31.0, ..., 3.808988764044944, 5.50561797752809,
        3.808988764044944],
       [-122.45, 37.77, 52.0, ..., 1.9859154929577465, 4.843505477308295,
        1.9859154929577465]], dtype=object)

五、数据预处理

# 获取训练数据（去除标签）
housing = strat_train_set.drop("median_house_value",axis=1)
# 获取训练数据的标签
housing_labels = strat_train_set['median_house_value'].copy()
# 获取测试集（去除标签）
test_dataset = strat_test_set.drop("median_house_value",axis=1)
# 获取训练集标签
test_labels = strat_test_set['median_house_value'].copy()

1. 处理缺失值

对于 total_bedrooms 属性的缺失，我们可以使用中位数进行填充

方法一

median_housing_value = housing['total_bedrooms'].median()
sample_incompletedata['total_bedrooms'].fillna(median_housing_value,inplace=True)

方法二

## Scikit-Learn 提供了容易的缺失值处理方式：SimpleImputer
from sklearn.impute import SimpleImputer

# 删除文字属性 ocean_proximity
housing_num = housing.drop('ocean_proximity',axis=1)

# 创建 imputer 对象
imputer = SimpleImputer(strategy="median")

# 为所有属性生成填充策略
imputer.fit(housing_num)

# 查看每个属性要替换的值
# imputer.statistics_

# 完成填充,结果是Numpy数组
X = imputer.transform(housing_num)

# 将Numpy数组转换回DataFrame格式
housing_tr = pd.DataFrame(X,columns=housing_num.columns)
housing_tr.info()

2. 处理文本和类别属性

LabelEncoder 转化器

# 将分类文本转换为数字形式
from sklearn.preprocessing import LabelEncoder
encoder = LabelEncoder()
housing_cateory = housing['ocean_proximity']
housing_cateory_encoded = encoder.fit_transform(housing_cateory)

# 查看转换器学习的映射关系，'<1H OCEAN' 是 0
# encoder.classes_
# 查看替换结果
# housing_cateory_encoded.reshape(-1,1)

3. one - hot 编码

OneHotEncoder 转换器

# 将数字形式的编码，转换为独热编码
from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder(categories='auto')
housing_cateory_1hot = encoder.fit_transform(housing_cateory_encoded.reshape(-1,1))
# 返回的结果是一个稀疏矩阵，以节省内存空间。可以调用.toarray()方法，将其转换为numpy数组。
housing_cateory_1hot.toarray()

array([[1., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 1.],
       ...,
       [0., 1., 0., 0., 0.],
       [1., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0.]])

LabelBinarizer 转换器

使用 LabelBinarizer 转换器，可以一次将文本转换为独热编码

from sklearn.preprocessing import LabelBinarizer
encoder = LabelBinarizer()
housing_cateor_1hot = encoder.fit_transform(housing_cateory)
housing_cateor_1hot

array([[1, 0, 0, 0, 0],
       [1, 0, 0, 0, 0],
       [0, 0, 0, 0, 1],
       ...,
       [0, 1, 0, 0, 0],
       [1, 0, 0, 0, 0],
       [0, 0, 0, 1, 0]])

4. 特征缩放

归一化与标准化

from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

num_pipline = Pipeline([
    ('imputer',SimpleImputer(strategy='median')),
    ('attribs_addr', CombineAttributesAdder()),
    ('std_scalar',StandardScaler()),
])

housing_num_tr = num_pipline.fit_transform(housing_num)
# 转化为 DataFrame 
columns = list(housing_num.columns) + ['rooms_per_household', 'population_per_household', 'bedrooms_per_room']
housing_num_tr = pd.DataFrame(housing_num_tr,columns=columns)
housing_num_tr.head()

5. 使用流水线处理

# 自定义转换器，将DataFrame对象转换为Numpy数组

from sklearn.base import BaseEstimator, TransformerMixin

class DataFrameSelector(BaseEstimator,TransformerMixin):
    def __init__(self,attribute_name):
        self.attribute_name = attribute_name
        
    def fit(self,X,Y=None):
        return self
    
    def transform(self,X):
        return X[self.attribute_name].values

# class V2toV1(BaseEstimator,TransformerMixin):
#     def fit(self,X,Y=None):
#         return self
    
#     def transform(self,X):
#         return X.reshape((-1,))

class MyLabelBinarizer(BaseEstimator,TransformerMixin):
    
    def __init__(self,encoder):
        self.encoder = encoder
        
    def fit(self,X,Y=None):
        self.encoder.fit(X)
        #返回self，确保在转换器中能够进行链式调用（例如调用transformer.fit(X).transform(X)）
        return self
    
    def transform(self,X, Y = None):
        return self.encoder.transform(X)
    

# FeatureUnion 类，可以将多个流水线合并在一起。

from sklearn.pipeline import FeatureUnion
num_attribs = list(housing_num.columns)
cat_attribs = ['ocean_proximity']
label_encoder = LabelBinarizer()


# 数值处理
num_pipline = Pipeline([
    ('selector',DataFrameSelector(num_attribs)),
    ('imputer',SimpleImputer(strategy='median')),
    ('attribs_addr',CombineAttributesAdder()),
    ('std_scaler',StandardScaler())
])

# 文本数据处理
cat_pipline = Pipeline([
    ('selector',DataFrameSelector(cat_attribs)),
    ('label_binarizer',MyLabelBinarizer(label_encoder))
])

full_pipelime = FeatureUnion(transformer_list=[
    ("num_pipeline", num_pipline),
    ("cat_pipeline", cat_pipline)
])

# 运行整条流水线
housing_prepared = full_pipelime.fit_transform(housing)
test_prepared = full_pipelime.fit_transform(test_dataset)

六、模型的建立

在进行模型的建立的时候，我们可以多尝试几种模型，最后可以选择一个效果比较好的，甚至可以组合一下他们。

1. 线性回归模型

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
# 训练模型
lin_reg.fit(housing_prepared,housing_labels)
y_predicted = lin_reg.predict(housing_prepared)

# 使用 mean_squared_error 来测量模型的 MSE
from sklearn.metrics import mean_squared_error

housing_predictions = lin_reg.predict(housing_prepared)
lin_mse = mean_squared_error(housing_labels, housing_predictions)

lin_rmse = np.sqrt(lin_mse)
lin_rmse

68628.19819848923

# 计算测试集是均方误差
test_predictions = lin_reg.predict(test_prepared)
test_lin_mse = mean_squared_error(test_labels, test_predictions)
test_lin_rmse = np.sqrt(test_lin_mse)
test_lin_rmse

66973.71087993948

2. 决策树模型

# 训练决策树模型
from sklearn.tree import DecisionTreeRegressor

tree_reg = DecisionTreeRegressor()
tree_reg.fit(housing_prepared,housing_labels)

# 计算决策树的均方误差
housing_predictions = tree_reg.predict(housing_prepared)
tree_mse = mean_squared_error(housing_labels,housing_predictions)
tree_rmse = np.sqrt(tree_mse)
tree_rmse

0.0

# 使用测试集数据计算决策树的均方误差
test_predictions = tree_reg.predict(test_prepared)
test_tree_mse = mean_squared_error(test_labels, test_predictions)
test_tree_rmse = np.sqrt(test_tree_mse)
test_tree_rmse

112797.81258965808

使用训练集数据计算的均方误差为 0 ，而使用测试集计算的均方误差却很大，这就是过拟合的问题，对于决策树来说，我们可以通过限制其树的深度,或者控制其每个节点最小分配数量来抑制过拟合的问题。

train_mse = []
test_mse = []
parameter_values = range(3,20)
for i in parameter_values:
    tree_reg = DecisionTreeRegressor(max_depth=i)
    tree_reg.fit(housing_prepared,housing_labels)
    housing_predictions = tree_reg.predict(housing_prepared)
    tree_mse = mean_squared_error(housing_labels,housing_predictions)
    tree_rmse = np.sqrt(tree_mse)
    train_mse.append(tree_rmse)
    test_predictions = tree_reg.predict(test_prepared)
    test_tree_mse = mean_squared_error(test_labels, test_predictions)
    test_tree_rmse = np.sqrt(test_tree_mse)   
    test_mse.append(test_tree_rmse)
    
plt.plot(parameter_values, train_mse, '-o', test_mse, '--o')
plt.show()

train_mse = []
test_mse = []
parameter_values = range(3,20)
for i in parameter_values:
    tree_reg = DecisionTreeRegressor(min_samples_leaf=i)
    tree_reg.fit(housing_prepared,housing_labels)
    housing_predictions = tree_reg.predict(housing_prepared)
    tree_mse = mean_squared_error(housing_labels,housing_predictions)
    tree_rmse = np.sqrt(tree_mse)
    train_mse.append(tree_rmse)
    test_predictions = tree_reg.predict(test_prepared)
    test_tree_mse = mean_squared_error(test_labels, test_predictions)
    test_tree_rmse = np.sqrt(test_tree_mse)   
    test_mse.append(test_tree_rmse)
    
plt.plot(parameter_values, train_mse, '-o', test_mse, '--o')
plt.show()

根据上述分析，我们最终选择参数：min_samples_leaf=15, max_depth=3，从下面的结果来看确实有一些提升

# 训练决策树模型
from sklearn.tree import DecisionTreeRegressor

tree_reg = DecisionTreeRegressor(min_samples_leaf=15, max_depth=3)
tree_reg.fit(housing_prepared,housing_labels)

# 计算决策树的均方误差
housing_predictions = tree_reg.predict(housing_prepared)
tree_mse = mean_squared_error(housing_labels,housing_predictions)
tree_rmse = np.sqrt(tree_mse)
tree_rmse

75116.29056699535

# 使用测试集数据计算决策树的均方误差
test_predictions = tree_reg.predict(test_prepared)
test_tree_mse = mean_squared_error(test_labels, test_predictions)
test_tree_rmse = np.sqrt(test_tree_mse)
test_tree_rmse

74841.17529820421

3. 随机森林

from sklearn.ensemble import RandomForestRegressor

forest_reg = RandomForestRegressor(n_estimators=20, max_depth=5)
forest_reg.fit(housing_prepared,housing_labels)

housing_predictions = tree_reg.predict(housing_prepared)
random_tree_mse = mean_squared_error(housing_labels,housing_predictions)
random_tree_rmse = np.sqrt(random_tree_mse)
random_tree_rmse

75116.29056699535

test_predictions = forest_reg.predict(test_prepared)
test_forest_mse = mean_squared_error(test_labels, test_predictions)
test_forest_rmse = np.sqrt(test_forest_mse)
test_forest_rmse

68024.4709392742

3. 模型泛化能力的对比

为了验证在不同数据下模型的表现效果，我们可以使用交叉验证的方法：

from sklearn.model_selection import cross_val_score

# Scikit-learn 的交叉验证倾向使用效用函数（值越大越好），而不是成本函数。所以这里使用负MSE。
scores = cross_val_score(tree_reg,housing_prepared,housing_labels, scoring='neg_mean_squared_error',cv=10)
rmse_scores = np.sqrt(-scores)

def display_scores(scores):
    print("Scores:",scores)
    print("Mean:", scores.mean())
    print("Std:",scores.std())

display_scores(rmse_scores)

## 查看线性模型的整体评分
lin_scores = cross_val_score(lin_reg, housing_prepared,housing_labels, scoring='neg_mean_squared_error',cv=10)
lin_rmse_scores = np.sqrt(-lin_scores)
lin_rmse_scores
display_scores(lin_rmse_scores)

forest_scores = cross_val_score(forest_reg, housing_prepared,housing_labels, scoring='neg_mean_squared_error',cv=10)
forest_rmse_scores = np.sqrt(-forest_scores)
forest_rmse_scores
display_scores(forest_rmse_scores)

从上面的三个模型的运行情况对比来看，准确度最高的是线性模型，交叉验证准确度最高的是决策树模型，但是这个模型不是很稳定，相比之下，随机森林泛化能力比较好。

七、总结

本文从最开始的数据获取，描述性统计分析、特征工程、再到最后的建模，建立了一个能估计加利福尼亚地区房屋价格的模型，本文较我之前关于数据分析的文章有一些新的突破，但是还是有很多地方有待改进。

1. 新的突破

• 使用了估计器、转换器、流水线，并根据自己的需要尝试实现了几个转换器，
• 在数据分析的过程中使用了 Tabeau 辅助分析

2. 有待改进

• 在特征工程部分做的很少，不知道如何进行
• 对模型的应用不是很熟练，缺少模型的误差来源分析与模型的改进，调参也不是很仔细，准确度还是太低