0. 写在前面
这是Boston房价数据分析的第二部分，中期的模型搭建与训练评估。大致的流程是，从原始论文的假设出发，先看看原作者对于该数据集的应用方式以及他所论证的问题。然后把文献的reference限制在相近的年份中，寻找在1978这个年代中使用该数据集的一些paper，这些paper的侧重点和研究方法会有不同，这里暂时先简单的做一下，因为这个项目还是需要机器学习的应用，如果纯写academic的话就太枯燥了，而且1978年的Boston，研究价值不大
为了以更加舒服的方式写这篇分析，可能会有中英文混用，等这个项目完成了有空可以详细翻译一下。其次是文中文献引用格式非常不标准，用markdown写有点麻烦，就随便引用了一下，保证可以找到原文即可。
1. 理论假设基础
1.1 原始数据集论文详解
这个数据集来自 Harrison and RubinfeldHarrison Jr, D., & Rubinfeld, D. L. (1978). Hedonic housing prices and the demand for clean air. 这篇论文的research goal是构建一个模型，用housing market data去measure the willingness to pay for clean air，所以本意上来说，数据集原始地应用场景是考量 housing market data和willingness to pay for clean air的关系，简单点说就是愿意花钱去买好房子=愿意花钱购买好的空气质量。
文章分为5个Sections (除过Introduction)
Section 2描述了一个four-step procedural model用来提出基本的实证研究。Step1是把空气质量作为享乐型房屋的属性，构建房屋价值模型; Step2是依据这个享乐型房屋价值模型，计算人们愿意为了 marginal change in air pollution 付出的意愿; Step3是通过先前的估计，来预估这种willingness在城市地区的模型，实际上就是愿意为了好的空气而付出的成本; Step4就是通过这个模型，来预估，当控制空气污染的策略生效，会给每单位房价带来的benefits (理解为，付出 $1来控制环境污染，那么购买房屋所花费的 $1会给房屋价值带来多少的提升)。
这个模型有一个前提假设，即房屋效用(U)是一系列属性的集合，而非只具有单一商品属性(housing is appropriately viewed as a bundle of attributes, rather than as a single commodity)，而air pollution是attributes中的一项 $a$。
$$
U(x,h)\
y = x + p(h) +T
$$
其中，$h$ 是房屋属性的集合即 ${h_1,h_2,…, h_n}$，$p(h)$是每个地方房屋价格的函数，$T$是运输成本，$x$是单位价格下复合私人产品的数量 (这里有些模糊，后面看看有没有具体解释)，$y$是budget。要使得 $U(x, h)$ 最大化，需要服从 $y$ 这个budget限制，即房屋价值不能超过其预算价值，而预算价值是由其属性构成的。
基于这个描述，Harrison提出了四点点关于享乐型住房的假设，
- 所有的顾客可以准确的感知不同地区房屋所具有的 $h$ 属性 (\向量，$\vec h$ )
- $ \vec h$ 的变化是连续的，所以 $ p(h) $ 房屋价格的变化是连续的，其一阶导和二阶导都是连续的
- 市场处于短期均衡
- 所有由于空间变化所产生的房屋属性 $ \vec h$ 的变化都可以被转换为房屋价格的变化
*这里需要注意，Harrison说明的很清楚，房屋属性与房屋价格的关系不一定是线性的 ( $p(h)$ relationship between housing attributes and house prices *need not to be linear)， 而非线性关系产生的原因是，市场不处于长期均衡 (Nonlinearities may exist in part because the market may not be in long-run equilibrium)，延申的去推断，根据假设3，这篇文章里使用的应该是线性关系。并且，房屋属性的集合不能被拆分重组为任意的属性集合 (housing attributes cannot to untied and repackaged to produce an arbitrary set of attributes at all locations) **
Step2是构建willingness to pay for marginal changes in air pollution, 一提到marginal，首先想到导数，所以在这部分，Harrison提出 “人们愿意为了小小地提升空气质量的付出 等同于 购买或租赁具有相同其他属性房屋所产生的投入的增加”，即 $[W_a(h)] = [p_a(h)]$
$$
W_a(h) = \frac {\partial U/ \partial(-a)} {\partial U/ \partial x} = \frac {\partial p(h)} {\partial (-a)} = p_a(h)
$$
Step3则是通过回归方法来估计 $W_a(h)$，这个step在section4。即当air pollution或者其他变量(household income)发生改变时，对于cleaner air的需求也会发生改变。在预估 $W_a(h)$ 之前需要先明确从可以获得的housing market data中，能不能识别出 $W_a(h)$ 。在这里，Harrison假设房屋的供求关系对于air pollution是锚定的，或者说sensitive是极低的，也就是说 $p(h)$ 受 $a$ 影响而非市场供给变化的影响。p.s 这种假设只是为了方便研究，市场对于房屋价格的影响在现实生活中是不可忽视的，在短期供需关系稳定的市场下，房价的波动可能也不会特别明显。thus, in my opinion, 一般情况下房价还是受到市场影响程度比较高(政策也是影响因素之一，在中国)，这与是否长期供求稳定关系不大，所以用纯线性关系来预测解释略有不妥。
Step4就是关于控制环境污染，与其带来的收益的关系讨论了，在section5有讨论。下面会具体结合paper与dataset来进行解读，对于论文中提到的分析方法着重研究。
首先，文章主要是在探讨空气质量和房屋价格，所以他的处理方式都是围绕着空气质量和房屋价格的关系来的。
NOX是衡量air pollution的重要指标，MEDV作为因变量的效果不如Log(MEDV)， 因为Log(MEDV)可以着重关注于NOX的非线性项。这里同时还有对NOX的转换问题，因为不同的计算方法会导致作用在房价的方式不同，论文采用的是 ${NOX}^2$ (原表达式是, ${NOX}^{p-1}/(p-1)$, p是NOX的参数) 。NOX呈现的负相关表现，是重点， 论文里还有一个PART也是和NOX一样，反映空气质量的。对于 ${NOX}^p$ 中 $p$ 值的选择，论文尝试了很多种，在原文Section5有详细的解释。
当删除DIS和RAD的时候，NOX的影响系数会从-0.0064上升到-0.0036，对此现象的原因是，Boston就业中心和公路附近的NOX比较高，删除这两项会削弱NOX的影响能力，结论是不应该删除。
当删除LSTAT时，NOX的影响系数会从-0.0064下降到-0.0081，对此现象的解释是，当删除LSTAT时，会把把NOX浓度归因于低层次人群导致的，而非其他复杂因素。
原文中的table1列举了5个等式，分为两类：Linear和Log-Log， 其中 $R^2$表现最好的是0.64，方程有两个：
$$
\log W = 1.05+0.78\log{Nox}+1.01\log {INC} - 0.24\log {PDU} \
\log W = 2.20+0.97\log {NOX} + 0.80\log{INC} \
- 0.03(Y_1)(\log{NOX})- 0.07(Y_2)(\log{NOX})
$$
而最终方程式的关系是：
$$
\log (MEDV) = a_1+a_2{RM^2}+a_3{AGE}+a_4\log {DIS} \
- a_5\log{RAD}+a_6{TAX}+a_7{PTRATIO}+a_8{(B-0.63)^2}\
- a_9\log{LSTAT}+ a_{10}{CRIM}+a_{11}{ZN}+a_{12}{INDUS}\
- a_{13}{CHAS}+ a_{14}{NOX}^p + \epsilon
$$
从这个方程中主要看到的是Harrison如何处理各个变量，使得最后的equation最佳，其中各个变量的含义和sklearn.dataset.load_boston()返回的数据集定义有些许不同，具体如下(原论文Table 4)
着重看一下与dataset表意不同的几个变量
- $B$ 在原数据集中等于 $1000(B_k-0.63)^2$, 其中 $B_k$才是黑人比例。 而在这里是 $B$ 就是黑人比例，处理在equation中是 $ (B-0.63)^2 $, 对比dataset的处理缩放了1000倍，这与第一部分中对于用 $B$ 反推 $B_k $ 并且缩放的判断基本相似。Harrison还解释了，为什么 $B$ 的形态是一个曲线：在 $B$ 的中低水平下，如果白人邻居对黑人视为不受欢迎的邻居类型， $B$ 发生增长将会带来housing price的负面影响。但是，市场歧视同时意味着当 $B$ 特别高的时候，housing price也会很高。 可以理解为，黑人在白人为主的社区中不受欢迎，但是在以黑人为主的社区中就不存在这样的问题。
LSTAT中需要注意，取对数的意义是：社会经济地位的区别在上层社会中的意义要大于下层阶级MEDV在论文中没有用 $’000 表示NOX在论文中的单位是 per 100 millionPART没有出现在 dataset中
在结束对原论文的分析后，可以发现 1. 如果要复现paper的研究，则不能删除变量，所以在传统的 linear model 中不讨论feature selection 2. 可以不用线性模型去做分析 3. paper的关注重点是air pollution 和 housing value，这种侧重决定了关于变量的处理是围绕 NOX来进行的，在后续模型中，如果采用不同的侧重点，则对features的处理方式也可以改变
1.2 其他论文的发现
主要是对于 $B$ 和 LSTAT的一些文章，没有详细的去看，引用的方式也很不标准。
Schafer, R. (1979). Racial discrimination in the Boston housing market. Journal of Urban Economics, 6(2), 176-196. 这篇文章提出：*黑人和低水平人群在房屋上的消费更高，白人愿意为了远离黑人社区而支付更高的购房价格，所以种族歧视现象在Boston Housing Market中存在的现象比较明显 *， 类似观点的还有Feins, J. D., & Bratt, R. G. (1983). Barred in Boston: Racial discrimination in housing. Journal of the American Planning Association, 49(3), 344-355.
基于这个观点，$B$ 可能不是像Harrison表述的是一个曲线，而应该是 $B_K$ 对 MEDV有正相关。所以在做模型的时候有不同的选择， 到底是用 $B$ 还是 $B_K$ 会对模型产生不同的影响。
2. 构建基本模型
终于到了要敲代码的地方了，在这里构建的模型主要会被分为两大类，第一，关于机器学习的模型；第二，关于神经网络的模型。在这里需要注意的是，机器学习的模型可以手撸，或者好好看看scikit-learn的文档，里面真的什么都有。
整个机器学习的模型设计思路比较简单，数据切割-输入-拟合-输出-真实值比对，其中，最重要的应该是对比真实值的结果，也就是拟合的程度。对比真实值，或者话说评估模型拟合程度的标准有很多，比如MSE均方误差，RMSE根均方误差， MAE平均绝对值误差，R²决定系数等等，在传统的统计学研究中，R²似乎是一个必须的选项。
2.1 机器学习模型
2.1.1 LiR线性回归
先来看看传统的未经处理的数据来进行LiR拟合
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import *
# 准备数据
data_boston = load_boston()
# 训练集和数据集准备 2/8分
x_train, x_test, y_train, y_test = train_test_split(data_boston.data, data_boston.target, train_size=0.8)
# 定义训练器
LiR = LinearRegression()
# 拟合
LiR.fit(x_train, y_train)
# 预测
y_predict = LiR.predict(x_test)
# 打分用的函数
LiR_standard_list = [mean_squared_error, mean_absolute_error, mean_gamma_deviance,
mean_poisson_deviance, mean_squared_log_error, mean_tweedie_deviance,
median_absolute_error, r2_score, max_error, explained_variance_score]
# 得分
def get_score(real, predict, standard_list, name):
for func in standard_list:
score = func(real, predict)
print(f"{name}模型下, {func.__name__} 的得分是：{score}")
if __name__ == '__main__':
get_score(y_test, y_predict, LiR_standard_list, "LiR")
"""
LiR模型下, mean_squared_error 的得分是：31.64095718595386
LiR模型下, mean_absolute_error 的得分是：3.7595557151678958
LiR模型下, mean_gamma_deviance 的得分是：0.10751142753841245
LiR模型下, mean_poisson_deviance 的得分是：1.471514913731351
LiR模型下, mean_squared_log_error 的得分是：0.0732933580733745
LiR模型下, mean_tweedie_deviance 的得分是：31.64095718595386
LiR模型下, median_absolute_error 的得分是：2.791386530941134
LiR模型下, r2_score 的得分是：0.6223503006075446
LiR模型下, max_error 的得分是：24.81682723887205
LiR模型下, explained_variance_score 的得分是：0.6284632981767782
"""在线性回归模型中，R²的得分很低，函数其实拟合的并不好，在自己测试的时候，通常要40000次左右的拟合，才会让score到90%以上，平均的得分只有70%左右，说明LiR这种模型多少还是有些不太给力。在原paper中也有提到用Log做处理，在这里模拟一下这个equation，流程是先处理输入值需要进行log的feature，然后把输出的y_predict进行log
"""
columns的名字和位置
RM 5 ,NOX 4 平方处理
DIS 7,RAD 8,LSTAT 12 Log处理
MEDV 0 先扩大1000，再log
B 11 缩小1000
"""
import numpy as np
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import *
# data准备
data_boston = load_boston()
# 预处理data，用新的数组存储
data = np.ndarray((506, 14), dtype=np.float64)
data[:, (4, 5)] = np.square(data_boston.data[:, (4, 5)])
data[:, (7, 8, 12)] = np.log(data_boston.data[:, (7, 8, 12)])
data[:, 13] = np.log(data_boston.target * 1000)
data[:, 11] = data_boston.data[:, 11] / 1000
data[:, (0, 1, 2, 3, 6, 9, 10)] = data_boston.data[:, (0, 1, 2, 3, 6, 9, 10)]
# 分割数据
x_train, x_test, y_train, y_test = train_test_split(data[:, 0:12], data[:, 13], train_size=0.8)
# 定义训练器
LiR = LinearRegression()
# 拟合
LiR.fit(x_train, y_train)
# 预测
y_predict = LiR.predict(x_test)
LiR_standard_list = [mean_squared_error, mean_absolute_error, mean_gamma_deviance,
mean_poisson_deviance, mean_squared_log_error, mean_tweedie_deviance,
median_absolute_error, r2_score, max_error,
explained_variance_score]
# 得分
def get_score(real, predict, standard_list, name):
for func in standard_list:
score = func(real, predict)
print(f"{name}模型下, {func.__name__} 的得分是：{score}")
if __name__ == '__main__':
get_score(y_test, y_predict, LiR_standard_list, "修正LiR")
"""
修正LiR模型下, mean_squared_error 的得分是：0.04546153727331139
修正LiR模型下, mean_absolute_error 的得分是：0.14697102771927917
修正LiR模型下, mean_gamma_deviance 的得分是：0.0005009416265240928
修正LiR模型下, mean_poisson_deviance 的得分是：0.004766181400402059
修正LiR模型下, mean_squared_log_error 的得分是：0.00041154977533747986
修正LiR模型下, mean_tweedie_deviance 的得分是：0.04546153727331139
修正LiR模型下, median_absolute_error 的得分是：0.11980174212924855
修正LiR模型下, r2_score 的得分是：0.7602699415038849
修正LiR模型下, max_error 的得分是：0.968815515553997
修正LiR模型下, explained_variance_score 的得分是：0.7625003112255664
"""可以发现，整个模型确实表现得更好了，但是这里也需要注意一下，整个get_score过程只有一次，是根据train_test_split所选取的x_train, y_train而来的，这里应该进行对折比对，循环切割n次数据集，得到总的平均得分，这样才能消除模型的偶然性positive或者negative得分，用平均得分来大致表示模型的拟合效果。
However，直接用最小二乘法会有个问题，那就是只有当矩阵 $X$ 是满秩的时候，才可以用最小二乘法。也就是多个因变量必须是相互独立的，如果相互之间关联较强，或者样本点比较少的时候，很可能造成$X$就不是满秩的，因为$X^{-1}$是不可逆的，而我们这个数据集中的共线性还是明显存在的，所以可以尝试使用带有L1,L2正则化的广义线性模型来处理，分别由LASSO(L1)和Ridge(L2)
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import *
from sklearn.metrics import *
data_boston = load_boston()
x_train, x_test, y_train, y_test = train_test_split(data_boston.data, data_boston.target, train_size=0.8)
lasso = LassoCV(cv=10, n_jobs=-1, n_alphas=2000)
ridge = RidgeCV(cv=10, alphas=(0.1, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0))
lasso.fit(x_train, y_train)
ridge.fit(x_train, y_train)
lasso_predict = lasso.predict(x_test)
ridge_predict = ridge.predict(x_test)
standard_lists = [mean_squared_error, mean_absolute_error, # mean_gamma_deviance,
mean_poisson_deviance, mean_squared_log_error, # mean_tweedie_deviance,
median_absolute_error, r2_score, max_error,
explained_variance_score]
# 得分
def get_score(real, predict, standard_list, name):
for func in standard_list:
score = func(real, predict)
print(f"{name}模型下, {func.__name__} 的得分是：{score}")
if __name__ == '__main__':
get_score(y_test, lasso_predict, standard_lists, "Lasso")
get_score(y_test, ridge_predict, standard_lists, "Ridge")
"""
注意，评分的时候需要把mean_gamma_deviance 和 mean_tweedie_deviance两个删除掉，不然会报错
Lasso模型下, mean_squared_error 的得分是：25.743721143836865
Lasso模型下, mean_absolute_error 的得分是：3.7328843155458724
Lasso模型下, mean_poisson_deviance 的得分是：0.9214281693677849
Lasso模型下, mean_squared_log_error 的得分是：0.03447203694086201
Lasso模型下, median_absolute_error 的得分是：2.831462338823245
Lasso模型下, r2_score 的得分是：0.6865925900801331
Lasso模型下, max_error 的得分是：17.894130269052802
Lasso模型下, explained_variance_score 的得分是：0.6870996145995208
Ridge模型下, mean_squared_error 的得分是：19.136565333389388
Ridge模型下, mean_absolute_error 的得分是：3.183047279909366
Ridge模型下, mean_poisson_deviance 的得分是：0.6738900372570746
Ridge模型下, mean_squared_log_error 的得分是：0.025707369207744177
Ridge模型下, median_absolute_error 的得分是：2.5099659672183936
Ridge模型下, r2_score 的得分是：0.7670289643680469
Ridge模型下, max_error 的得分是：15.629044054083195
Ridge模型下, explained_variance_score 的得分是：0.7685006048843454
"""
详细介绍一下各个score的算法
mean_squared_error: $\text{MSE}(y, \hat{y}) = \frac{1}{n_\text{samples}} \sum_{i=0}^{n_\text{samples} - 1} (y_i - \hat{y}_i)^2.$ Best
mean_absolute_error: $\text{MAE}(y, \hat{y}) = \frac{1}{n_{\text{samples}}} \sum_{i=0}^{n_{\text{samples}}-1} \left| y_i - \hat{y}_i \right|.$ Best
mean_tweedie_deviance; mean_poisson_deviance (power=1) , mean_gamma_deviance (power=2), mean_squared_error (power=0):
$$
\begin{split}\text{D}(y, \hat{y}) = \frac{1}{n_\text{samples}}
\sum_{i=0}^{n_\text{samples} - 1}
\begin{cases}
(y_i-\hat{y}i)^2, & \text{for }p=0\text{ (Normal)}\
2(y_i \log(y/\hat{y}_i) + \hat{y}_i - y_i), & \text{for}p=1\text{ (Poisson)}\
2(\log(\hat{y}_i/y_i) + y_i/\hat{y}_i - 1), & \text{for}p=2\text{ (Gamma)}\
2\left(\frac{\max(y_i,0)^{2-p}}{(1-p)(2-p)}-
\frac{y,\hat{y}^{1-p}_i}{1-p}+\frac{\hat{y}^{2-p}_i}{2-p}\right),
& \text{otherwise}
\end{cases}\end{split}
$$mean_squared_log_error: $\text{MSLE}(y, \hat{y}) = \frac{1}{n_\text{samples}} \sum{i=0}^{n_\text{samples} - 1} (\log_e (1 + y_i) - \log_e (1 + \hat{y}_i) )^2.$
median_absolute_error: $\text{MedAE}(y, \hat{y}) = \text{median}(\mid y_1 - \hat{y}_1 \mid, \ldots, \mid y_n - \hat{y}_n \mid).$
R²: $R^2(y, \hat{y}) = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}i)^2}{\sum{i=1}^{n} (y_i - \bar{y})^2}$
max_error : $\text{Max Error}(y, \hat{y}) = max(| y_i - \hat{y}_i |)$
在参考的博文中，发现LiR模型和kNN模型的得分差距不大，但是在SVR和LiR对比时，LiR普遍在R²得分较高，而SVR在MAE等得分较高，所以就把SVR作为第二个机器学习的模型来进行尝试。
2.1.2 SVR支持向量机回归
传统回归方法当且仅当回归f(x)完全等于y时才认为预测正确，如线性回归中常用$(f(x)−y)^2$来计算其损失。而支持向量回归则认为只要$f(x)$与$y$偏离程度不要太大，既可以认为预测正确，不用计算损失，具体的，就是设置阈值α,只计算 $|f(x)−y|>α$ 的数据点的loss
简单点说，LiR下预测出来的是一条线，误差是散点与线的之间的差，从数学角度来说：点无大小，线无粗细，所以平面上的散点完全落在直线上的概率很低，这也是LiR不能完全准确的预测出$y$值的原因，也是导致过拟合或者欠拟合的一种可能。在SVR的视角下，根据kernel的不同，预测的是一个区域，凡是在区域内的点都“可以”被视为无误差的，加之C惩罚因子的存在，既能考虑离群点，又能防止过拟合，很棒棒哦。
from sklearn.svm import SVR
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.metrics import *
data_boston = load_boston()
x_train, x_test, y_train, y_test = train_test_split(data_boston.data, data_boston.target, train_size=0.8)
poly_svr = SVR(kernel='poly')
liner_svr = SVR(kernel='linear')
rbf_svr = SVR(kernel='rbf')
svr_list = [poly_svr, liner_svr, rbf_svr]
standard_list = [mean_squared_error, mean_absolute_error, median_absolute_error,
v_measure_score, r2_score, max_error, explained_variance_score]
def get_score(real, predict, standard_list, name):
for sc in standard_list:
score = sc(real, predict)
print(f"{name}模型下, {sc.__name__} 的得分是：{score}")
if __name__ == '__main__':
for func in svr_list:
func.fit(x_train, y_train)
y_predict = func.predict(x_test)
get_score(y_test, y_predict, standard_list, func.kernel)
"""
poly模型下, mean_squared_error 的得分是：37.50059531965942
poly模型下, mean_absolute_error 的得分是：4.069654995770053
poly模型下, median_absolute_error 的得分是：2.703214692506716
poly模型下, v_measure_score 的得分是：0.9654092985785466
poly模型下, r2_score 的得分是：0.25681850585358934
poly模型下, max_error 的得分是：35.04311418113312
poly模型下, explained_variance_score 的得分是：0.25991748287715377
linear模型下, mean_squared_error 的得分是：17.365140241461017
linear模型下, mean_absolute_error 的得分是：2.564467175077125
linear模型下, median_absolute_error 的得分是：1.905132283343118
linear模型下, v_measure_score 的得分是：0.9654092985785466
linear模型下, r2_score 的得分是：0.6558601067342158
linear模型下, max_error 的得分是：28.458002879364713
linear模型下, explained_variance_score 的得分是：0.656437977890237
rbf模型下, mean_squared_error 的得分是：38.158738286894895
rbf模型下, mean_absolute_error 的得分是：4.197313783839424
rbf模型下, median_absolute_error 的得分是：3.359745270874889
rbf模型下, v_measure_score 的得分是：0.9654092985785466
rbf模型下, r2_score 的得分是：0.24377552161340066
rbf模型下, max_error 的得分是：34.603823044964656
rbf模型下, explained_variance_score 的得分是：0.2452859318065379
"""
可以看到，我自己做出来的数据和参考博文的情况不太一致，一般情况下来说，rbf的kernel应该是拟合情况最好的，可是MAE普遍较大，而且R²真的是表现不佳。
2.1.3 GBM 梯度提升机
GBM是一种Boosting的实例算法，其核心思想是分阶段集合弱预测模型以产生新的模型, 这种算法梯度增强以迭代方式将弱的“学习器”组合为一个强学习器。最简单的解释是在最小二乘回归中，通过最小化均方误差(MSE)，“教”模型预测实数值 $ \hat{y}$。再简单一点说，线性回归的一个基本假设是其残差之和为0，即残差应该在0左右随机分布，GBM的思想就是迭代优化那些离散程度大的点，使其拟合程度更好，而这个过程也就是在优化MSE使其最小。
而同家族Boosting算法下，还有GBDT基于回归树的梯度提升算法，也叫提升树算法BRT，其核心原理和GBM类似，是迭代多棵回归树来共同决策。当采用MSE时，每一棵回归树学习的是之前所有树的结论和残差，拟合得到一个当前的残差回归树，提升树即是整个迭代过程生成的回归树的累加。
在先前的博客中看到有人使用LightGBM(GBDT升级版)，并且有不错的结果，那么这里就看一下LGBM的表现。
from lightgbm.sklearn import LGBMRegressor
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.metrics import *
data_boston = load_boston()
x_train, x_test, y_train, y_test = train_test_split(data_boston.data, data_boston.target, train_size=0.8)
LGBM = LGBMRegressor(objective="regression",
boosting_type='gbdt',
num_iterations=300,
learning_rate=0.01,
num_leaves=31,
num_threads=8,
n_estimators=30)
LGBM.fit(x_train, y_train)
y_predict = LGBM.predict(x_test)
standard_list = [mean_squared_error, mean_absolute_error, mean_gamma_deviance,
mean_poisson_deviance, mean_squared_log_error, mean_tweedie_deviance,
median_absolute_error, v_measure_score, r2_score, max_error,
explained_variance_score]
# 得分
def get_score(real, predict, standard_list, name):
for func in standard_list:
score = func(real, predict)
print(f"{name}模型下, {func.__name__} 的得分是：{score}")
if __name__ == '__main__':
get_score(y_test, y_predict, standard_list, "LGBM")
"""
LGBM模型下, mean_squared_error 的得分是：9.664227163518742
LGBM模型下, mean_absolute_error 的得分是：2.083257201677075
LGBM模型下, mean_gamma_deviance 的得分是：0.027396120698409937
LGBM模型下, mean_poisson_deviance 的得分是：0.48653665987789796
LGBM模型下, mean_squared_log_error 的得分是：0.02481375736532234
LGBM模型下, mean_tweedie_deviance 的得分是：9.664227163518742
LGBM模型下, median_absolute_error 的得分是：1.3402386997906603
LGBM模型下, v_measure_score 的得分是：0.9799501732117103
LGBM模型下, r2_score 的得分是：0.858467809515831
LGBM模型下, max_error 的得分是：12.839740252866987
LGBM模型下, explained_variance_score 的得分是：0.8586600270975745
"""整体下来LGBM的得分确实很好，这可能也是为什么比赛里很多人都使用GBM算法的原因了。
2.2 神经网络模型
原本计划参考用NN和BP网络去做回归，可是发现模型没有解释意义，构建一个单输入和单输出，最后不用激活函数的神经网络，虽然可以应用于连续回归问题，但是其解释意义太弱了，不如单纯的机器学习模型。所以暂时不用网络去做了，等自己再学习学习，看看有没有好的解释办法。
3. 写在最后
这部分东西拖了很久才做，中途好多事情影响。1. 写了个链家的爬虫，结果IP池做的不好，爬取成功率只有80%左右，然后又去优化了IP池，改成了高匿池。2. 把原来不用的笔记本翻出来，做了个arch的虚拟机，把爬虫放上去调试了一段时间 3. 受疫情影响，求学就业之路不顺，自闭了很久。 3. 机器学习很多算法层面的东西学的不扎实，补了很多课，还是觉得不够好，又开始自闭了 4. hexo疯狂报错，真的想自己写个框架了。。。。。。
慢慢来慢慢来
