端到端的机器学习项目（求加州住房价格中位数)

下载数据

当数据会定期变化时，写一个脚本，在需要的时候直接运行（或者设置一个定期自动运行的计划任务），或者需要在多台机器上安装数据集，自动化函数非常好用

注：因为网站的问题，所以没有科学上网的话可能不能执行
所以数据集放在这里
加州房价数据集

import os
import tarfile
import urllib

DOWNLOAD_ROOT = "http://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):
    os.makedirs(housing_path, exist_ok=True)
    tgz_path = os.path.join(housing_path, "housing.tgz")
    # print(tgz_path)
    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()  # 只有第一次运行时需要运行这个代码，之后将这个代码注释掉即可

使用pandas加载数据

import pandas as pd
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.head()

运行结果：

查看数据

housing.info()  # 从中可以看出ocean_proximity为非数值数据，total_bedrooms是缺值数据

运行结果：
从上述运行结果可以看到不同指标的Non-Null Count值不完全相同，在这里可以看出total_bedrooms对应的值要小，说明他存在空值，并且可以看到Dtype（数据类型）中除去ocean_proximity是object之外均是float64，为数值类型。

housing["ocean_proximity"].value_counts()  # 非数值部分的数据

运行结果：

housing.describe()  # 只统计数值数据

%matplotlib inline  # 注意，这行代码只适用于jupyter notebook
import matplotlib.pyplot as plt
housing.hist(bins=55, figsize=(20, 15))
plt.show()

数据数值分布

数据划分

1.随机划分

import numpy as np

def split_train_test(data, test_ratio):
    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(housing, 0.2) 

优点：便于理解
缺点：这种方法虽然可以进行划分。但是每次的运行结果都不一样
解决方法：（1）第一次运行之后就保存下来（当新数据加入后需要重新进行）（2）设置随机种子如np.random.seed(42)（3）哈希

2.哈希方法

from zlib import crc32

def test_set_check(identifier, test_ratio):
    return crc32(np.int64(identifier)) & 0xffffffff < test_ratio * 2 ** 32  # 哈希函数

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

housing_with_id = housing.reset_index()  # 添加索引, housing_with_id是添加之后的数据
train_set, test_set = split_train_test_by_id(housing_with_id, 0.2, "index")

注：使用行索引作为标识符必须确保只在末尾添加数据，并且不会删除任何行，在本问题中因为经纬度不会重复，所以可以利用经纬度来进行唯一标识。

housing_with_id["id"] = housing["longitude"] * 1000 + housing["latitude"]
train_set, test_set = split_train_test_by_id(housing_with_id, 0.2, "id")

3.scikit-Learn划分子集

from sklearn.model_selection import train_test_split

train_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)

下面来思考一个问题，如果说数据分成多个层次，每个层次数据量不同，如果按照随机划分的话数据量大的数据明显比数据量小的数据有优势。所以要采用分层
4.纯随机的抽样会出现显著偏差，需要分层取数据

设定标签并分层

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]) 
# 分层不宜过多，这里分了5层

housing.head()  # income_cat存的是pd.cut之后的标签

运行结果（由于数据过长，图只截了一部分，具体结果请大家自己运行一下）：

housing["income_cat"].hist()  # 显示各个标签的个数

按照收入类别划分训练集和测试集

from sklearn.model_selection import StratifiedShuffleSplit

split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)
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]

strat_train_set["income_cat"].value_counts() / len(strat_train_set)

运行结果

strat_test_set["income_cat"].value_counts() / len(strat_test_set)

运行结果：

去除income_cat

for set_ in (strat_train_set, strat_test_set):
    set_.drop("income_cat", axis=1, inplace=True)
strat_train_set.describe() # 可以看到这里没有了income_cat，数据恢复了原样

通过数据可视化观测数据

创建副本（为了防止原数据被修改）

housing = strat_train_set.copy() # 现在来讲能探索的只是训练集，为了不损害数据集，先要创建一个副本

地理位置可视化

housing.plot(kind="scatter", x="longitude", y="latitude")  # 用经度当x坐标，纬度为y坐标，画散点图

运行结果：

看起来和加州的地图很像（因为图本来就是按照经纬度画的），但是在这个图中除了得到房子位置信息之外我们得不到其他信息了。假如我们要突出高密度区域，我们可以在函数中加入alpha参数（参数大小具体视情况而定），这样可以突出高密度地区，弱化低密度地区（抓大放小）。

housing.plot(kind="scatter", x="longitude", y="latitude", alpha=0.1)  # 突出高密度区域的更好的可视化

housing.plot(kind="scatter", x="longitude", y="latitude", alpha=0.4,
             s=housing["population"]/100, label="population",figsize=(10, 7),  # s表示的是圆点的大小，这里是使用人口数量作为圆的大小
             c="median_house_value", cmap=plt.get_cmap("jet"), colorbar=True,  # 颜色代表价格（选项c）, jet定义颜色表（选项cmap）来可视化
             # 这里是用颜色来表示平均房价
             # 颜色范围从蓝（低）到红（高）
            )
plt.legend()

# 图像表明房价和地理位置和人口密度相关，通常使用聚类算法检测主集群，然后再为每个集群中心添加一个新的衡量临近距离的特征。
# 海洋临近度是个很有用的特征，但是在加州沿海地区的房价也不是很高，所以这个简单的规则也不是万能的

运行结果：

现在对于数据已经有了一个基本的了解了，接下来开始寻找数据的相关性（在机器学习中，更偏好使用和预测目标相关性较高的指标）

# 由于数据量不是很大，可以使用corr()的方式计算每对属性之间的标准相关系数（也称皮尔逊r）
corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)  # 数据靠近0说明两者之间没有线性关系

运行结果：

相关系数越接近1相关性越强，可以看到这里边除去平均房价本身是1，最高的是平均收入。

另一种查看相关性的方式

from pandas.plotting import scatter_matrix

attributes = ["median_house_value", "median_income", "total_rooms", "housing_median_age"]
scatter_matrix(housing[attributes], figsize=(12, 8))

运行结果：

从上面两种查看相关系数的方法可以看到平均收入最接近线性，我们把他分离出来单独看一下

housing.plot(kind="scatter", x="median_income", y="median_house_value", alpha=0.1)

运行结果：

可以看到出现了一些直线，如50万，45万，35万，28万…为了防止将来算法出现这种怪异的数据，可能要删除一部分数据，那一部分可能是异常值或者噪声。

试验不同属性的组合

因为有些原始数据相关性较小，而经过处理后相关性会增加即

# 在准备给机器学习算法输入数据前，最后一件事是尝试各种属性的组合
housing["rooms_per_houshold"] = housing["total_rooms"] / housing["households"]
housing["bedrooms_per_room"] = housing["total_bedrooms"] / housing["total_rooms"]
housing["popluation_per_houshold"] = housing["population"] / housing["households"]
corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)

运行结果：

可以看到bedrooms_per_room比rooms_per_household和population_per_househoud与房价中位数的相关性更好

数据准备

准备数据时应该编写函数，而不是手动进行
1.可以在任意数据集上重现这些转换
2.可以逐渐建立转换函数的函数库，可以在以后的项目中使用
3.可以在实时系统中使用这些函数来转换新数据，再输入给出算法
4.可以尝试是、多种转换方式，查看哪种最佳

现在我们先回到一个干净的训练集

housing = strat_train_set.drop("median_house_value", axis=1)
housing_labels = strat_train_set["median_house_value"].copy()

housing_labels.head()

运行结果：

数据清洗

大部分的机器学习算法都无法在缺失的特征上工作，如total_bedrooms属性有部分值的缺失

方法：
1.放弃这些区域的数据
2.放弃整个属性
3.将缺失的值设置为某个值（0， 平均数或中位数）

实现方式：
方法1：housing.dropna(subset=[“total_bedrooms”])
方法2：housing.drop(“total_bedrooms”, axis=1)
方法3：median = housing[“total_bedrooms”].median()
housing[“total_bedrooms”].fillna(median, inplace=True)

scikit-learn提供了一个非常容易上手的类处理缺失值

from sklearn.impute import SimpleImputer

# 新建一个SimppleImputer实例，指定要用属性的中位数替换该属性的缺失值
imputer = SimpleImputer(strategy="median")

# 由于中位数值只能在数值属性上计算，所以我们需要创建一个没有文本属性ocean_proximity数据副本
housing_num = housing.drop("ocean_proximity", axis=1)

# 使用fit()方法，将imputer实例适配到训练数据
imputer.fit(housing_num)
# 这里imputer仅仅只是计算了每个属性的中位数值，并将结果存储在其实例变量statistics_中，
# 虽然只有total_bedrooms这个属性存在缺失，但是我们没办法确定系统启动后新数据是否一定不存在任何缺值，
# 稳妥起见将imputer应用于所有的数值属性值

imputer.statistics_

运行结果：

housing_num.median().values

运行结果：

用imputer将缺失值替换成中位数值从而完成训练集转换

X = imputer.transform(housing_num)  # 这是一个numpy格式的数据

housing_tr = pd.DataFrame(X, columns=housing_num.columns, index=housing_num.index)  # 因为X得出的是numpy格式的，要将其转换为DataFrame

housing_tr.head()  # 与之前的housing对比发现少了ocean_proximity（因为imputer只能用于数值的计算，在之前的处理中已经将其去掉了）

处理文本和分类属性

# 到目前为止只处理数值属性，现在看一下文本
housing_cat = housing[["ocean_proximity"]]
housing_cat.head(10)  # 默认参数是5，现在选用的参数是10，所以有10条

运行结果：

直接使用OrdinalEncoder将文本转换成数值

# 他是有限可能的取值，现在由文本转换为数字
from sklearn.preprocessing import OrdinalEncoder

ordinal_encoder = OrdinalEncoder()
housing_cat_encoded = ordinal_encoder.fit_transform(housing_cat)
housing_cat_encoded[:10]

ordinal_encoder.categories_

运行结果：

转换成one_hot编码

# 转换为one_hot编码
from sklearn.preprocessing import OneHotEncoder

cat_encoder = OneHotEncoder()
housing_cat_1hot = cat_encoder.fit_transform(housing_cat)  # 输出的是一个scipy的稀疏矩阵而不是numpy数组（只存储非0）,这样可以防止存储过多的0元

转换成array格式

housing_cat_1hot.toarray()

运行结果：

cat_encoder.categories_

运行结果：

自定义转换器

创建一个类，应用fit()(返回self)、transform()、fit_transform()

from sklearn.base import BaseEstimator, TransformerMixin

rooms_ix, bedrooms_ix, population_ix, households_ix = 3, 4, 5, 6

class CombinedAttributesAdder(BaseEstimator, TransformerMixin):
    # 添加TransformerMixin作为基类可直接得到fit_transform
    # 添加BaseEstimator作为基类（并在构造函数中避免*args和**kargs方法），
    # 可以额外获得两种非常有用的自动调整超参数的方法（get_params()和set_params）
    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):
        rooms_per_household = X[:, rooms_ix] / X[:, households_ix]
        population_per_household = X[:, population_ix] / X[:, households_ix]
        if self.add_bedrooms_per_room:
            bedroom_per_room = X[:, bedrooms_ix] / X[:,rooms_ix]
            return np.c_[X, rooms_per_household, population_per_household, bedroom_per_room]
        else:
            return np.c_[X, rooms_per_household, population_per_household]

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

转换流水线

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

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

housing_num_tr = num_pipline.fit_transform(housing_num)

转换后结果

训练和评估训练集

线性回归

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[:5]
some_data_prepared = full_pipeline.transform(some_data)
print("Predictions:", lin_reg.predict(some_data_prepared))
print("Labels:", list(some_labels))

运行结果：

均方差误差

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)
print(lin_rmse)

决策树回归

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_predictions, housing_labels)
tree_rmse = np.sqrt(tree_mse)
tree_rmse

运行结果：
0.0

使用交叉验证更好地进行评估

1.划分训练集和验证集
2.K折

交叉验证决策树

from sklearn.model_selection import cross_val_score

scores = cross_val_score(tree_reg, housing_prepared, housing_labels, scoring="neg_mean_squared_error", cv=10)
tree_rmse_score = np.sqrt(-scores)
tree_rmse_score

def display_scores(scores):
    print("Scores:", scores)
    print("Mean:", scores.mean())
    print("Standard deviation:", scores.std())
display_scores(tree_rmse_score)

交叉验证线性回归

lin_scores = cross_val_score(lin_reg, housing_prepared, housing_labels, scoring="neg_mean_squared_error", cv=10)
lin_rmse = np.sqrt(-lin_scores)
display_scores(lin_rmse)

交叉验证随机森林

from sklearn.ensemble import RandomForestRegressor

forest_reg = RandomForestRegressor()
forest_reg.fit(housing_prepared, housing_labels)
housing_predictions = forest_reg.predict(housing_prepared)
forest_mse = mean_squared_error(housing_predictions, housing_labels)
forest_rmse = np.sqrt(forest_mse)
print(forest_rmse)

运行结果：
18746.865963012133

在这里插入代码片

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)

模型保存和加载

import joblib
joblib.dump(lin_reg, "line_reg.pkl")
my_model_loaded = joblib.load("line_reg.pkl")

模型微调

网格搜索

from sklearn.model_selection import GridSearchCV

# 调整随机森林的超参数
param_grid = [
    {'n_estimators':[3, 10, 30], 'max_features':[2, 4, 6, 8]},
    {'bootstrap':[False], 'n_estimators':[3, 10], 'max_features':[2, 3, 4]}
]
forest_reg = RandomForestRegressor()
grid_search = GridSearchCV(forest_reg, param_grid, scoring='neg_mean_squared_error', return_train_score=True)
grid_search.fit(housing_prepared, housing_labels)

运行结果：

grid_search.best_params_

grid_search.best_estimator_

cvres = grid_search.cv_results_
cvres.keys()

for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
    print(np.sqrt(-mean_score), params)

随机搜索：超参数范围较大时优先使用RandomizeSearchCV

分析最佳模型及误差

# 随机森林中重要程度
feature_importances = grid_search.best_estimator_.feature_importances_
feature_importances

extra_attribs = ["room_per_hhold", "pop_per_hhold", "bedrooms_per_room"]
cat_encoder = full_pipeline.named_transformers_["cat"]
cat_one_hot_attribs = list(cat_encoder.categories_[0])
print(num_attribs)
print(extra_attribs)
print(cat_one_hot_attribs)
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)
print(final_mse)
print(final_rmse)

运行结果：

from scipy import stats

confidence = 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)))