1.模型目标

预测某一区域的房价中位数

2.选择框架

有监督学习任务：训练集中的每个实例都有标签(该区域的房价中位数)
回归任务：因为你要对某个值进行预测。更具体地说，这是一个多重回归问题，因为系统要使用多个特征进行预测（使用区域的人口、收入中位数等）。这也是一元回归问题，因为我们仅尝试预测每个区域的单个值。
简单的批量学习应该就能胜任：我们没有一个连续的数据流不断流进系统，所以不需要针对变化的数据做出特别调整。批量学习：当系统投入生产环境后，学习就停止。在线学习：系统发布了，还需要根据输入的训练数据逐步累积学习成果。

3.选择性能指标

均方根误差（RMSE）(Root Mean Square Error) --- 回归问题的首选性能指标

m表示数据集中的实例数
$x^{(i)}$ 是数据集中第i个实例的所有特征值的向量， $y^{(i)}$ 是其标签（该实例的期望输出值）
h是系统的预测函数，也称为假设。当给系统输入一个实例的特征向量 $x^{(i)}$ 时，它会为该实例输出一个预测值 $y\widehat{}^{(i)}$ =h( $x^{(i)}$ )
X是一个矩阵，包含数据集中所有实例的所有特征值。每个实例占一行，第 i 行等于 $x^{(i)}$ 的转置
RMSE（X，h）表示使用预测函数h在一组示例中测量的成本函数

平均绝对误差（MAE）(Mean Absolute Error)

我们将包含n个元素的向量v的 $l_{k}$ 范数定义为 $\left \| v \right \|_{k}=(|v_{0}|^{k}+|v_{1}|^{k}+...+|v_{n}|^{k})^{1/k}$

RMSE对应 $l_{2}$ 范数，又称欧几里得范数，其实就是几何空间中两点的距离
MAE对应 $l_{1}$ 范数，又称曼哈顿范数。如果你只能沿着正交的城市街区行动，它能够测量城市中两点之间的距离。
范数指标越高，它越关注大值而忽略小值。这就是RMSE对异常值比MAE更敏感的原因。

4.获取数据

import os
import tarfile
import urllib.request
import pandas as pdDOWNLOAD_ROOT = "https://raw.githubusercontent.com/ageron/handson-ml2/master/"
HOUSING_PATH = os.path.join("datasets", "housing")
HOUSING_URL = DOWNLOAD_ROOT + "datasets/housing/housing.tgz"def fetch_housing_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()def load_housing_data(housing_path=HOUSING_PATH):csv_path = os.path.join(housing_path, "housing.csv")return pd.read_csv(csv_path)fetch_housing_data() # 下载房价数据
housing = load_housing_data() # 获取数据对象,返回一个包含所有数据的pandas DataFrame对象

5.分析数据

①返回dataframe最上面的五行，可以发现共有十列

housing.head()

②快速获取数据集的简单描述，特别是总行数、每个属性的类型和非空值的数量

housing.info()

数据集中包含20640个实例
total_bedrooms这个属性只有20433个非空值，这意味着有207个区域缺失这个特征

所有属性的字段都是数字，除了ocean_proximity，它的类型是object，查看一下有多少种分类存在

 housing["ocean_proximity"].value_counts()

③通过describe()方法可以显示数值属性的描述

std行显示的是标准差（用来测量数值的离散程度）
关于25%、50%和75%，以housing_median_age为例，它的含义是：25%的值小于18，50%的值小于29，75%的值小于37

④绘制每个数值属性的直方图

%matplotlib inline # 指定Jupyter作为Matplotlib的后端
import matplotlib.pyplot as plt
# x轴数据集中房价的值，纵轴表示这个值出现的次数
# bins表示直方图的“柱”的个数,每个“柱”的值=该“柱”跨越的所有x对应的y值之和
housing.hist(bins=50, figsize=(20,15))
plt.show()

收入中位数为万美元
许多直方图都表现出重尾分布(Heavy-tailed distribution)：图形在中位数右侧的延伸比左侧要远得多。这可能会导致某些机器学习算法难以检测模式。
稍后我们会尝试一些转化方法，将这些属性转化为更偏向钟形的分布

6.创建测试集

随机选择一些实例，通常是数据集的20%

纯随机抽样

from sklearn.model_selection import train_test_splittrain_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)

分层抽样
随机抽样很可能带来误差，假设我们要对一个地区的人进行抽样调查，该地区的男女比例是3:1，那么抽取的样本应该保持这一比例。
同理，如果我们被告知：预测房价中位数，收入中位数是一个非常重要的属性
我们先看看收入中位数的直方图：大多数的值聚集在1.5～6左右

创建5个收入类别属性income_cat（用1～5来做标签），0～1.5是类别1，1.5～3是类别2，以此类推

import numpy as np
housing["income_cat"] = pd.cut(housing["median_income"],bins=[0., 1.5, 3.0, 4.5, 6., np.inf], labels=[1, 2, 3, 4, 5])
housing["income_cat"].hist()

根据收入类别进行分层抽样

from sklearn.model_selection import StratifiedShuffleSplitsplit = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42) # n_splits表示生成一组train/test对
# split(housing,housing["income_cat"])表示按照housing["income_cat"]各种值的比例在housing中抽取样本,样本的形式为数组下标
for train_index, test_index in split.split(housing, housing["income_cat"]):strat_train_set = housing.loc[train_index]strat_test_set = housing.loc[test_index]

检验：查看原数据集、随机test集、分层test集根据income_cat分类情况，并计算误差

def income_cat_proportions(data):return data["income_cat"].value_counts() / len(data)train_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)compare_props = pd.DataFrame({"Overall": income_cat_proportions(housing),"Stratified": income_cat_proportions(strat_test_set),"Random": income_cat_proportions(test_set),
}).sort_index()
compare_props["Rand. %error"] = 100 * compare_props["Random"] / compare_props["Overall"] - 100
compare_props["Strat. %error"] = 100 * compare_props["Stratified"] / compare_props["Overall"] - 100compare_props

删除income_cat属性

for set_ in (strat_train_set, strat_test_set):set_.drop("income_cat", axis=1, inplace=True)

7.可视化数据探索

首先，把测试集放在一边，你能探索的只有训练集。如果训练集非常庞大，你可以抽样一个探索集。

explore_set = strat_train_set.copy()

①由于存在地理位置信息（经度和纬度），因此建立一个各区域的分布图来看看

# 通过调整alpha，可以更清楚地看出高密度数据点的位置
explore_set.plot(kind="scatter", x="longitude", y="latitude", alpha=0.1)

②现在，再来看看房价

# 每个圆的半径大小代表了每个区域的人口数量（选项s）
# 颜色代表价格（选项c），颜色范围从蓝（低）到红（高）
explore_set.plot(kind="scatter", x="longitude", y="latitude", alpha=0.4,s=explore_set["population"]/100, label="population", figsize=(10,7),c="median_house_value", cmap=plt.get_cmap("jet"), colorbar=True,sharex=False)

③寻找相关性

可以使用corr()方法计算每对属性之间的标准相关系数（也称为皮尔逊相关系数）

# 看看每个属性与房价中位数的相关性分别是多少
# 相关系数的范围从-1变化到1。越接近1，表示有越强的正相关
# 系数靠近0则说明二者之间没有线性相关性
corr_matrix = explore_set.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)

由此可见，收入中位数median_income确实和房价关系很大

④试验不同属性的组合

如果你不知道一个区域有多少个家庭，单看total_rooms(房间总数)、total_bedrooms(卧室总数)似乎没啥意义，所以我们可以考虑组合新的属性

explore_set["rooms_per_household"] = explore_set["total_rooms"]/housing["households"] # 平均每户人家拥有的房间数目
explore_set["bedrooms_per_room"] = explore_set["total_bedrooms"]/housing["total_rooms"] # 房间总数和卧室总数的比例
explore_set["population_per_household"]=explore_set["population"]/housing["households"] # 平均每户人家有几口人
# 查看相关矩阵
corr_matrix = explore_set.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)

8.数据准备

先准备一份训练数据和对应的标签

train_data = strat_train_set.drop("median_house_value", axis=1) # 去掉"median_house_value"这一列作为训练集
train_labels = strat_train_set["median_house_value"].copy() # "median_house_value

①数据清理

前面我们已经注意到total_bedrooms属性有部分值缺失，所以我们要解决它。我们有三种选择：
1.放弃缺失属性所在的整个实例
2.放弃整个属性
3.将缺失的值设置为某个值（0、平均数或者中位数等）

train_data.dropna(subset=["total_bedrooms"])    # option 1
train_data.drop("total_bedrooms", axis=1)       # option 2
median = train_data["total_bedrooms"].median()  # option 3
train_data["total_bedrooms"].fillna(median, inplace=True)

Scikit-Learn提供了一个SimpleImputer类来处理缺失值。

from sklearn.impute import SimpleImputer
# 创建一个SimpleImputer实例，指定你要用属性的中位数值替换该属性的缺失值
imputer = SimpleImputer(strategy="median")
# 由于中位数值只能在数值属性上计算，所以我们需要创建一个没有文本属性ocean_proximity的数据副本
train_data_num = train_data.drop("ocean_proximity", axis=1)
# 使用fit()方法将imputer实例适配到训练数据
# imputer仅仅只是计算了每个属性的中位数值，并将结果存储在其实例变量statistics_中
# 也就是说，imputer.statistics_和train_data_num.median().values的值是一样的
imputer.fit(train_data_num)
# 利用imputer将缺失值替换成中位数值
X = imputer.transform(train_data_num)
# 结果是一个NumPy数组,转换为pandasDataFrame
housing_tr = pd.DataFrame(X, columns=train_data_num.columns, index=train_data_num.index)

②处理文本和分类属性

在此数据集中，ocean_proximity是文本属性，但它不是任意文本，它的取值是有限的，每个值代表一个类别，因此，此属性是分类属性。

train_data_cat = train_data[["ocean_proximity"]]
train_data_cat.head(10)

使用Scikit-Learn的OrdinalEncoder类，将这些类别从文本转到数字。

from sklearn.preprocessing import OrdinalEncoderordinal_encoder = OrdinalEncoder()
train_data_cat_encoded = ordinal_encoder.fit_transform(train_data_cat)
train_data_cat_encoded[:10]
# 可以使用Categories_实例变量获取类别列表
ordinal_encoder.categories_

整数编码的一个问题是，机器学习算法会认为两个相近的值比两个离得较远的值更为相似一些。
在某些情况下这是对的（对一些有序类别，像“坏”“平均”“好”“优秀”）
但是，对ocean_proximity而言并非如此（例如，类别0和类别4就比类别0和类别1的相似度更高）
解决：one-hot编码

from sklearn.preprocessing import OneHotEncodercat_encoder = OneHotEncoder()
# 这里返回的是一个SciPy稀疏矩阵，该稀疏矩阵仅存储非零元素,因为0太多了,占内存
train_data_cat_1hot = cat_encoder.fit_transform(train_data_cat)
# 你也可以将它转为numpy矩阵
train_data_cat_1hot.toarray()

当有很多这样的类别属性需要转换，one-hot编码会导致输入特征很多，会减慢训练并降低性能。
所以，你也可以考虑用相关的数字特征代替类别输入。
例如，你可以用与海洋的距离来替换ocean_proximity特征。

虽然Scikit-Learn提供了许多有用的转换器，但是你仍然需要为一些诸如自定义清理操作或组合特定属性等任务编写自己的转换器。这里限于篇幅就不展开了。

9.选择和训练模型

数据预处理：对上面的数据操作做一个汇总

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.compose import ColumnTransformer# column index
rooms_ix, bedrooms_ix, population_ix, households_ix = 3, 4, 5, 6# 自定义一个属性组合的转换器
class CombinedAttributesAdder(BaseEstimator, TransformerMixin):def __init__(self, add_bedrooms_per_room=True): # no *args or **kargsself.add_bedrooms_per_room = add_bedrooms_per_roomdef fit(self, X, y=None):return self  # 什么也不做def transform(self, X):rooms_per_household = X[:, rooms_ix] / X[:, households_ix]population_per_household = X[:, population_ix] / X[:, households_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]# Pipeline函数创建了一个管道，用于集成一些数据预处理的操作
num_pipeline = Pipeline([('imputer', SimpleImputer(strategy="median")), # 将缺失值替换为中位数('attribs_adder', CombinedAttributesAdder()),('std_scaler', StandardScaler()), # 数据归一化和标准化])num_attribs = list(train_data_num) # 所有数字属性的属性名称
cat_attribs = ["ocean_proximity"] # 分类属性名称# 自定义对所有数据处理的管道
full_pipeline = ColumnTransformer([("num", num_pipeline, num_attribs),("cat", OneHotEncoder(), cat_attribs), # 分类属性进行one-hot编码])housing_prepared = full_pipeline.fit_transform(train_data) # train_data表示分层抽样后的训练数据集strat_train_set去掉median_house_value一列后的数据集
housing_labels = train_labels # 标签

线性回归模型

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_errorlin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels) # 训练一个线性回归模型# 算一下均方根误差
housing_predictions = lin_reg.predict(housing_prepared)
lin_mse = mean_squared_error(housing_labels, housing_predictions)
lin_rmse = np.sqrt(lin_mse)
lin_rmse # 68627.87390018745

显然，RMSE的值有点大，说明该模型对训练数据欠拟合，想要修正欠拟合，可以通过选择更强大的模型，或为算法训练提供更好的特征。

决策树

from sklearn.tree import DecisionTreeRegressortree_reg = DecisionTreeRegressor(random_state=42)
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

以上，说明该模型在训练集上表现良好，但有可能模型过拟合了，我们可以用K-折交叉验证功能来验证一下(模型没搞好，建议先不触碰测试集)
以10-折交叉验证为例：将训练集随机分割成10个不同的子集，每个子集称为一个折叠，然后对决策树模型进行10次训练和评估——每次挑选1个折叠进行评估，使用另外的9个折叠进行训练。
产生的结果是一个包含10次评估分数的数组

from sklearn.model_selection import cross_val_scorescores = cross_val_score(tree_reg, housing_prepared, housing_labels, scoring="neg_mean_squared_error", cv=10)
tree_rmse_scores = np.sqrt(-scores)def display_scores(scores):print("Scores:", scores)print("Mean:", scores.mean())print("Standard deviation:", scores.std()) # 精确度(标准差)display_scores(tree_rmse_scores)
# 输出（显然过拟合了）:
Scores: [72831.45749112 69973.18438322 69528.56551415 72517.78229792 69145.50006909 79094.74123727 68960.045444 73344.50225684 69826.02473916 71077.09753998]
Mean: 71629.89009727491
Standard deviation: 2914.035468468928

随机森林

from sklearn.ensemble import RandomForestRegressorforest_reg = RandomForestRegressor(n_estimators=100, random_state=42)
forest_reg.fit(housing_prepared, housing_labels)
# 算一下均方根误差
housing_predictions = forest_reg.predict(housing_prepared)
forest_mse = mean_squared_error(housing_labels, housing_predictions)
forest_rmse = np.sqrt(forest_mse)
forest_rmse # 18650.698705770003

# 交叉验证
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)
display_scores(forest_rmse_scores)
# 输出
Scores: [51559.63379638 48737.57100062 47210.51269766 51875.21247297 47577.50470123 51863.27467888 52746.34645573 50065.1762751 48664.66818196 54055.90894609]
Mean: 50435.58092066179
Standard deviation: 2203.3381412764606

但是，训练集上的分数仍然远低于验证集，这意味着该模型仍然对训练集过拟合。
过拟合的可能解决方案：简化模型、约束模型（即使其正规化），或获得更多的训练数据。

10.微调模型

可以手动调整超参数，也可以借助Scikit-Learn的GridSearchCV（网格搜索）

from sklearn.model_selection import GridSearchCVparam_grid = [# 首先评估这3×4=12种超参数组合{'n_estimators': [3, 10, 30], 'max_features': [2, 4, 6, 8]},# 接着，尝试这2×3=6种超参数组合{'bootstrap': [False], 'n_estimators': [3, 10], 'max_features': [2, 3, 4]},]forest_reg = RandomForestRegressor(random_state=42)
# 对每个模型进行5次训练（因为使用的是5-折交叉验证），也就是说，总共会完成18×5=90次训练
grid_search = GridSearchCV(forest_reg, param_grid, cv=5, scoring='neg_mean_squared_error', return_train_score=True)
grid_search.fit(housing_prepared, housing_labels)

grid_search.best_params_ # 找到最佳的参数组合

输出：{'max_features': 8, 'n_estimators': 30}

grid_search.best_estimator_ # 直接得到最好的估算器
# 输出评估分数
cvres = grid_search.cv_results_
for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):print(np.sqrt(-mean_score), params)

综上，最佳解决方案：将超参数max_features设置为8，将超参数n_estimators设置为30。
这个组合的RMSE分数为49 898，略优于之前的模型。

随机搜索

如果探索的组合数量较少，那么网格搜索是一种不错的方法。
但是当超参数的搜索范围较大时，通常会优先选择使用RandomizedSearchCV。
这个类用起来与GridSearchCV类大致相同，但它不会尝试所有可能的组合，而是在每次迭代中为每个超参数选择一个随机值，然后对一定数量的随机组合进行评估。

集成方法

将表现最优的模型组合起来。组合（或“集成”）方法通常比最佳的单一模型更好（就像随机森林比其所依赖的任何单个决策树模型更好一样），特别是当单一模型会产生不同类型误差时更是如此。

分析最佳模型及其误差

在进行准确预测时，RandomForestRegressor可以指出每个属性的相对重要程度
将这些重要性分数显示在对应的属性名称旁边

feature_importances = grid_search.best_estimator_.feature_importances_
extra_attribs = ["rooms_per_hhold", "pop_per_hhold", "bedrooms_per_room"]
cat_encoder = full_pipeline.named_transformers_["cat"]
cat_one_hot_attribs = list(cat_encoder.categories_[0])
attributes = num_attribs + extra_attribs + cat_one_hot_attribs
sorted(zip(feature_importances, attributes), reverse=True)

有了这些信息，你可以尝试删除一些不太有用的特征

11.通过测试集评估系统

final_model = grid_search.best_estimator_X_test = strat_test_set.drop("median_house_value", axis=1)
y_test = strat_test_set["median_house_value"].copy()X_test_prepared = full_pipeline.transform(X_test) # 预处理后的数据
final_predictions = final_model.predict(X_test_prepared)final_mse = mean_squared_error(y_test, final_predictions)
final_rmse = np.sqrt(final_mse)
final_rmse # 47873.26095812988

在某些情况下，泛化误差(rmse)的这种点估计不足以说服你启动生产环境，你需要知道这个估计的精确度。
为此，你可以使用scipy.stats.t.interval()计算泛化误差的95%置信区间

from scipy import statsconfidence = 0.95
squared_errors = (final_predictions - y_test) ** 2
np.sqrt(stats.t.interval(confidence, len(squared_errors) - 1, loc=squared_errors.mean(), scale=stats.sem(squared_errors)))

输出：array([45893.36082829, 49774.46796717])