《蜥蜴书》_讲义及源码解读_02

完整的机器学习项目

使用开源数据集（真实的）

https://archive.ics.uci.edu/ml/ UC Irvine Machine Learning Repository

https://www.kaggle.com/datasets Kaggle datasets

https://registry.opendata.aws/ aws.amazon.com/cn/s3/ 需要注册

https://tianchi.aliyun.com/dataset

https://www.datatang.com/dataset/all

https://github.com/awesomedata/awesome-public-datasets

https://github.com/wangqingbaidu/Dr.Sure

http://dataportals.org/

放两个入门的：

http://yann.lecun.com/exdb/mnist/

http://www.vision.caltech.edu/Image_Datasets/Caltech101/

示例项目，依据加州房产价格数据集（基于 1990 年加州普查的数据），预测房价。

划定问题

你的模型的输出（预测一个街区的房价中位数）会传入后面的系统。这一整套系统可以确定某个街区进行投资分析和投资决策，直接影响利润。很重要。

如图，一系列的数据处理组件被称为数据流水线。流水线在机器学习系统中很常见，因为有许多数据要处理和转换。

调研目前的状况：

• 专家队伍（人工）处理数据，费时费钱，结果不理想，误差大概有15%，确定BaseLine!
• 这是一个典型的监督学习任务，因为你要使用的是有标签的训练样本（每个实例都有预定的产出，即街区的房价中位数）。
• 这是一个典型的回归任务，因为你要预测一个值。同时是一个多变量回归问题，因为系统要使用多个变量进行预测（要使用街区的人口，收入中位数等等）
• 没有连续的数据流进入系统，没有特别需求需要对数据变动作出快速适应。数据量不大可以放到内存中，因此批量学习就够了。

选择性能指标

回归问题的典型指标是：均方根误差（RMSE）。均方根误差测量的是系统预测误差的标准差。
$$RMSE(X,h)=\sqrt{\frac{1}{m}\sum _{i=1}^m(h(x^i)-y^i{)}^2}$$$$m:数据集中实例样本数量$$$$x^i:第i个实例样本的所有特征值向量$$$$y^i:第i个实例样本的标签。监督学习$$$$X:数据集中所有实例样本的所有特征值的矩阵（不包含标签）$$$$h：是预测函数（也成为假设）,输入一个样本实例的特征向量，输出一个预测值,即$$$${\hat{y}}^i=h(x^i)$$

在有些情况下（实例样本数据中有许多异常时），可能需要使用：平均绝对误差（MAE）
$$MAE(X,h)=\frac{1}{m}\sum _{i=1}^m|h(x^i)-y^i|$$

RMSE 和MAE 都是测量预测值和目标值两个向量距离的方法。有多种测量距离的方法，或范数：

RMSE对应的是欧几里得范数（欧式距离），也成为L2范数，均方根误差，记为：$||\cdot |{|}_2$

MAE对应的是曼哈顿范数，也称为L1范数，平均绝对误差，记为：$||\cdot |{|}_1$

范数指标越高，它越关注大值而忽略小值，这就是RMSE对异常值比MAE更敏感的原因

检查假设

和流水线上下游沟通，进一步核实前期的假设是否成立，比如：确定下游是不是需要：预测的价格，如果下游需要的是房产的分类呢，比如：“昂贵、中端、廉价”，这个问题就变成分类问题，而不是回归问题了。

准备开发环境并下载数据

编写函数：实现下载，解压的自动化。

import os
import tarfile   #压缩文件处理
import urllib.request   #网络请求
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

DOWNLOAD_ROOT = "https://raw.githubusercontent.com/ageron/handson-ml2/master/"
HOUSING_PATH = os.path.join("datasets", "housing")
HOUSING_URL = DOWNLOAD_ROOT + "datasets/housing/housing.tgz"

#PATH文件夹  ：当前目录下的：datasets/housing   
#HOUSING_URL：
#https://raw.githubusercontent.com/ageron/handson-ml2/master/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() # 关闭压缩文件

    
fetch_housing_data()    #下载、解压数据文件

#读取文件为PD的dataframe格式，已知数据文件名为：housing.csv
def load_housing_data(housing_path=HOUSING_PATH):
    csv_path = os.path.join(housing_path, "housing.csv")
    return pd.read_csv(csv_path)

housing = load_housing_data()  #调用函数读取文件housing.csv
print(housing.shape)  #显示数据形状  (20640, 10)
housing.head()  #　显示开始的5行数据 ，正：东经,北纬　　负：西经,南纬

查看数据结构

housing.info()  #  显示数据结构信息  

RangeIndex: 20640 entries, 0 to 20639
Data columns (total 10 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   longitude           20640 non-null  float64  经度
 1   latitude            20640 non-null  float64  纬度
 2   housing_median_age  20640 non-null  float64  房龄中位数
 3   total_rooms         20640 non-null  float64  总房间数
 4   total_bedrooms      20433 non-null  float64  总卧室数
 5   population          20640 non-null  float64  人口数量
 6   households          20640 non-null  float64  家庭数（户数）
 7   median_income       20640 non-null  float64  收入中位数
 8   median_house_value  20640 non-null  float64  房价中位数
 9   ocean_proximity     20640 non-null  object   邻近海洋
dtypes: float64(9), object(1)

查看所有数值型特征的统计指标

housing.describe()  # 显示数据的整体统计分析情况

查看非数值特征邻近海洋 的分类汇总

housing["ocean_proximity"].value_counts()
#显示数据中"ocean_proximity"特征的分类计数结果

<1H OCEAN     9136
INLAND        6551
NEAR OCEAN    2658
NEAR BAY      2290
ISLAND           5
Name: ocean_proximity, dtype: int64

绘制直方图，查看数据分布情况

housing.hist(bins=50, figsize=(20,15))  # 50个样本实例为1组
plt.show()

创建训练集、测试集

分割数据集，创建训练集和测试集的比例一般为8：2

洗牌1：普通的随机抽样方式

from sklearn.model_selection import train_test_split

train_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)
# 
#X_train,X_test, y_train, y_test =train_test_split(X,y,test_size, random_state)

洗牌2：本例根据median_income 的分布，将收入中位数特征分成5类，根据类别进行分层抽样。

from sklearn.model_selection import StratifiedShuffleSplit

#增加income_cat 特征，依据前面直方图观察到的分布分5类，标记为1-5
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])

ss = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)
#n_splits=1  默认10，这里只产生1组洗牌效果，分别对应训练集索引和测试机的索引

#提取n组分层抽样结果，对于本例其实只有1组，
for train_index, test_index in ss.split(housing, housing["income_cat"]):
    strat_train_set = housing.loc[train_index]
    strat_test_set = housing.loc[test_index]

#检查效果，测试集income_cat的占比比例和原始数据集一致。
strat_test_set["income_cat"].value_counts() / len(strat_test_set)  
housing["income_cat"].value_counts() / len(housing)

#当然，拿到训练集和测试集后，可以删掉临时的income_cat特征。参考df.drop()函数
for set_ in (strat_train_set, strat_test_set):
    set_.drop("income_cat", axis=1, inplace=True)

数据探索和可视化

地理数据角度可视化

#拿到分层抽样的训练集的一个副本，以免损伤训练集：
housing = strat_train_set.copy()
#因为存在地理信息（纬度和经度），创建一个所有街区的散点图来数据可视化
housing.plot(kind="scatter", x="longitude", y="latitude")

#将alpha设为 0.1，透明度，可以更容易看出数据点的密度
housing.plot(kind="scatter", x="longitude", y="latitude", alpha=0.1)

以房价角度可视化

每个圈的半径表示街区的人口（选项s），颜色代表价格（选项c）。

我们用预先定义的名为jet的颜色图（选项cmap），它的范围是从蓝色（低价）到红色（高价）

housing.plot(kind="scatter", x="longitude", y="latitude", alpha=0.4,
             s=housing["population"]/100, label="population", figsize=(10,7),
             c="median_house_value", cmap=plt.get_cmap("jet"), colorbar=True,
             sharex=False)
plt.legend()

#对本例，叠加上加州地图，来查看房价分布情况，更直观，尤其看到近邻海岸对房价的影响。

寻找相关性

因为数据集并不是非常大，你可以很容易地使用corr()方法计算出每对属性间的标准相关系数（standard correlation coefficient，也称作皮尔逊相关系数）

import seaborn as sns
from pandas.plotting import scatter_matrix
#拿到协方差矩阵，相关系数矩阵。参考pd.corr()函数
corr_matrix = housing.corr()

sns.heatmap(corr_matrix)  #流行的seaborn库的热图

#根据相关系数矩阵，查看相关性，拿到相关度相对大的4个特征，绘制散点矩阵
attributes = ["median_house_value", "median_income", "total_rooms",
              "housing_median_age"]



scatter_matrix(housing[attributes], figsize=(12, 8))


#两种可视化相关系数矩阵，结论：和房价相关性最大的是收入

关系数的范围是 -1 到 1。当接近 1 时，意味强正相关；当相关系数接近 -1 时，意味强负相关；相关系数接近 0，意味没有线性相关性。图 2-14 展示了相关系数在横轴和纵轴之间的不同图形。

实验不同属性的组合

给算法准备数据之前，你需要做的最后一件事是尝试多种属性组合

对于本例：

• 如果你不知道某个街区有多少户，该街区的总房间数就没什么用。
• 你真正需要的是每户有几个房间。 总房间数/户数
• 相似的，总卧室数也不重要：你可能需要将其与房间数进行比较。卧室数/总房间数
• 每户的人口数也是一个有趣的属性组合，人口数/户数。让我们来创建这些新的属性（特征）：

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()
sns.heatmap(corr_matrix)
corr_matrix["median_house_value"].sort_values(ascending=False)
#可看出，组合属性rooms_per_household 相关性，比他们未组合前各自的相关性要高很多。

机器学习算法的数据准备

分割下训练数据和标签

housing = strat_train_set.drop("median_house_value", axis=1) 
# drop labels for training set  drop(),不会破坏原数据 ，inplace默认为false
housing_labels = strat_train_set["median_house_value"].copy()

数据清洗

机器学习算法不能处理缺失的特征，因此先创建一些函数来处理特征缺失的问题。

本例中，属性total_bedrooms有一些缺失值。有三个解决选项：

• 放弃对应的街区（整行）；

• 放弃整个属性（整列）；

• 进行赋值（0、平均值、中位数等等）

• Scikit-Learn提供了一个方便的类来处理缺失值：Imputer。首先，需要创建一个Imputer实例，指定用某属性的中位数来替换该属性所有的缺失值

• from sklearn.impute import SimpleImputer
  #拿到SimpleImputer的实例，策略是"median"
  imputer = SimpleImputer(strategy="median") 
  
  housing_num = housing.drop("ocean_proximity", axis=1)
  #先去掉ocean_proximity，因为中位数只能用于数值型数据的计算
  imputer.fit(housing_num)  
  imputer.statistics_  #输出结果数组，是每特征的中位数
  housing_num.median().values  #检查一下是否正确
  
  X = imputer.transform(housing_num)  #先fit，后transform ,结果是个numpy数组
  housing_tr = pd.DataFrame(X, columns=housing_num.columns,index=housing.index)
  #将numpy数组转换成pd,  现在得到的就是所有数值型特征中，用中位数填补缺失值的训练数据集。
  
  

处理文本和分类特征

机器学习算法中，需要文本标签转换为数值。本例中，ocean_proximity 特征是文本属性。

方法一：直接编码

housing_cat = housing[["ocean_proximity"]]
housing_cat.head(10)

from sklearn.preprocessing import OrdinalEncoder

ordinal_encoder = OrdinalEncoder()
housing_cat_encoded = ordinal_encoder.fit_transform(housing_cat)
housing_cat_encoded[:10]
#用fit_transform，直接编码
#array([[0.],
#       [0.],
#       [4.],
#       [1.],
#       [0.],
#       [1.],
#       [0.],
#       [1.],
#       [0.],
#       [0.]])
ordinal_encoder.categories_

#[array(['<1H OCEAN', 'INLAND', 'ISLAND', 'NEAR BAY', 'NEAR OCEAN'],
#       dtype=object)]

方法2：使用onehot 编码

from sklearn.preprocessing import OneHotEncoder
cat_encoder = OneHotEncoder(sparse=False)
housing_cat_1hot = cat_encoder.fit_transform(housing_cat)
housing_cat_1hot
housing_cat_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.]])

cat_encoder.categories_

特征缩放

两种常见的方法可以让所有的属性（特征）有相同的度量

• 归一化(Min-Max scaling): 通过减去最小值，然后再除以最大值与最小值的差值，取值范围[0，1]
• 标准化(Standardization): 首先减去平均值，然后除以方差。标准化受到异常值的影响很小。

综上，转换流水线

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

num_pipeline = Pipeline([
        ('imputer', SimpleImputer(strategy="median")),
        ('attribs_adder', CombinedAttributesAdder()),
        ('std_scaler', StandardScaler()),
    ])

housing_num_tr = num_pipeline.fit_transform(housing_num)
#中位数填补空缺，自定义转换器（增加属性），特征缩放， 依次流水，得到housing_num_tr

#新版sklearn 提供了新的这：ColumnTransformer
from sklearn.compose import ColumnTransformer

num_attribs = list(housing_num)
cat_attribs = ["ocean_proximity"]

full_pipeline = ColumnTransformer([
        ("num", num_pipeline, num_attribs),
        ("cat", OneHotEncoder(), cat_attribs),
    ])
housing_prepared = full_pipeline.fit_transform(housing)

选择并训练模型

在训练集上训练和评估

使用普通的线性回归

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels) #训练！


some_data = housing.iloc[:5]
some_labels = housing_labels.iloc[:5]
some_data_prepared = full_pipeline.transform(some_data) 

#预测并观察对比前5个实例，这里还是训练集，不碰测试集。
print("Predictions:", lin_reg.predict(some_data_prepared))
print("Labels:", list(some_labels))

#Predictions: [210644.6 317768.8 210956.4  59218.9 189747.5]
#Labels: [286600.0, 340600.0, 196900.0, 46300.0, 254500.0]

#计算下rmse  
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)  # 68628.198  ,欠拟合。

#计算下mae
from sklearn.metrics import mean_absolute_error

lin_mae = mean_absolute_error(housing_labels, housing_predictions)
# 49439.896  ,欠拟合！

使用决策树模型试试

from sklearn.tree import DecisionTreeRegressor

tree_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)
# 0，均方误差为0？ ，神奇的过拟合！

使用交叉验证更好的对决策树模型进行评估

from sklearn.model_selection import cross_val_score
#K-折交叉验证，CV=10,训练集分成10个子集，用9个进行训练，1个进行评估
scores = 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)
#得到10个评估结果:
Scores: [70194.33 66855.16 72432.58 70758.73 ... 71231.65]
Mean: 71407.68766037929
Standard deviation: 2439.4345041191004
#决策树的评分大约是 71407.687，通常波动有 ±2439 
               

结论：比普通线性回归的结果还要糟糕，真实的反映。

再来试试随机森林模型

from sklearn.ensemble import RandomForestRegressor

forest_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)
#18603.515021376355  用训练集上的数据预测后计算mse


from sklearn.model_selection import cross_val_score

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)
#交叉验证后的mse
#Mean: 50182.303100336096
#Standard deviation: 2097.0810550985693

结论，比前几模型好，但还是过拟合严重。交叉验证的mse 远远大于训练集的mse

机器学习的十八般武艺，还可以在搞一搞SVM 支撑向量机模型，等等。

微调模型

假设挑选了一个靠谱的模型，或一个模型列表。可以手工调整超参数，直到找到一个好的超参数组合，但低效、鼓噪，耗时到生无可恋。

网格搜索

from sklearn.model_selection import GridSearchCV
#网格参数：列表数组
param_grid = [
    # 尝试 12 (3×4) 种超参数组合
    {'n_estimators': [3, 10, 30], 'max_features': [2, 4, 6, 8]},
    # 然后尝试 6 (2×3) 种bootstrap 为fasle时的超参数组合
    {'bootstrap': [False], 'n_estimators': [3, 10], 'max_features': [2, 3, 4]},
  ]
#选择前面表现较好的随机森林模型
forest_reg = RandomForestRegressor(random_state=42)
# CV=5,
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}

#当然还可以修正param_grid网格参数，继续搜索更优的超参数组合

随机搜索

当搜索次数相对较少时，如上面$（3\times 4+2\times 3）\ast 5=90$次，适合使用网格搜索。

但是当超参数的搜索空间很大时，最好使用RandomizedSearchCV。

它不是尝试所有可能的组合，而是通过选择每个超参数的一个随机值的特定数量的随机组合。

这个方法有两个优点：

• 如果你让随机搜索运行，比如 1000 次，它会探索每个超参数的 1000 个不同的值（而不是像网格搜索那样，只搜索每个超参数的几个值）。
• 你可以方便地通过设定搜索次数，控制超参数搜索的计算量。

from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint

param_distribs = {
        'n_estimators': randint(low=1, high=200),
        'max_features': randint(low=1, high=8),
    }

forest_reg = RandomForestRegressor(random_state=42)
rnd_search = RandomizedSearchCV(forest_reg, param_distributions=param_distribs,
                                n_iter=10, cv=5, scoring='neg_mean_squared_error', 		random_state=42)
rnd_search.fit(housing_prepared, housing_labels)

#可以同样拿到相对最好的超参数组合
#49280.9449827171 {'max_features': 7, 'n_estimators': 122}

分析最佳模型及其误差

通过网格搜索还可以拿到特征重要性的指标，已确定特征对模型的重要性。

feature_importances = grid_search.best_estimator_.feature_importances_
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)

通过测试集评估系统

最后，碰一碰测试集，使用测试集数据验证下模型的优劣。

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)

启动、监控、维护系统

保存模型

方法一：pickle模块 import pickle

方法二：sklearn内部的joblib from sklearn.externals import joblib

发布模型

法一：上传到网站，随web服务器启动时加载，通过客户端访问请求出发服务端调用预测方法predict()

法二：上传到云，或者说部署到云。

监控与维护模型

编写监控代码，以固定间隔检测系统的实时表现，当发生性能下降时触发报警

引入人工分析者对模型采样评估，分析者可以是领域专家，或是众包平台的工人

定期用新数据训练模型。你应该尽可能自动化这个过程