Sklearn-秘籍-二-Sklearn 秘籍（二）第二章处理线性模型本章包括以下主题：线性回归模型评估线性回

Sklearn 秘籍（二）

协议：CC BY-NC-SA 4.0

第二章处理线性模型

作者：Trent Hauck

译者：muxuezi

协议：CC BY-NC-SA 4.0

本章包括以下主题：

简介

线性模型是统计学和机器学习的基础。很多方法都利用变量的线性组合描述数据之间的关系。通常都要花费很大精力做各种变换，目的就是为了让数据可以描述成一种线性组合形式。

本章，我们将从最简单的数据直线拟合模型到分类模型，最后介绍贝叶斯岭回归。

2.1 线性回归模型

现在，我们来做一些建模！我们从最简单的线性回归（Linear regression）开始。线性回归是最早的也是最基本的模型——把数据拟合成一条直线。

Getting ready

boston数据集很适合用来演示线性回归。boston数据集包含了波士顿地区的房屋价格中位数。还有一些可能会影响房价的因素，比如犯罪率（crime rate）。

首先，让我们加载数据：

from sklearn import datasets
boston = datasets.load_boston()

How to do it...

实际上，用scikit-learn的线性回归非常简单，其API和前面介绍的模型一样。

首先，导入LinearRegression类创建一个对象：

from sklearn.linear_model import LinearRegression
lr = LinearRegression()

现在，再把自变量和因变量传给LinearRegression的fit方法：

lr.fit(boston.data, boston.target)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

2.2 评估线性回归模型

在这个主题中，我们将介绍回归模型拟合数据的效果。上一个主题我们拟合了数据，但是并没太关注拟合的效果。每当拟合工作做完之后，我们应该问的第一个问题就是“拟合的效果如何？”本主题将回答这个问题。

Getting ready

我们还用上一主题里的lr对象和boston数据集。lr对象已经拟合过数据，现在有许多方法可以用。

from sklearn import datasets
boston = datasets.load_boston()
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(boston.data, boston.target)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

predictions = lr.predict(boston.data)

How to do it...

我们可以看到一些简单的量度（metris）和图形。让我们看看上一章的残差图：

%matplotlib inline
from matplotlib import pyplot as plt
f, ax = plt.subplots(figsize=(7, 5))
f.tight_layout()
ax.hist(boston.target - predictions,bins=40, label='Residuals Linear', color='b', alpha=.5);
ax.set_title("Histogram of Residuals")
ax.legend(loc='best');

png

如果你用IPython Notebook，就用%matplotlib inline命令在网页中显示matplotlib图形。如果你不用，就用f.savefig('myfig.png')保存图形，以备使用。

画图的库是matplotlib，并非本书重点，但是可视化效果非常好。

和之前介绍的一样，误差项服从均值为0的正态分布。残差就是误差，所以这个图也应噶近似正态分布。看起来拟合挺好的，只是有点偏。我们计算一下残差的均值，应该很接近0：

import numpy as np
np.mean(boston.target - predictions)

6.0382090193051989e-16

另一个值得看的图是Q-Q图（分位数概率分布），我们用Scipy来实现图形，因为它内置这个概率分布图的方法：

from scipy.stats import probplot
f = plt.figure(figsize=(7, 5))
ax = f.add_subplot(111)
probplot(boston.target - predictions, plot=ax);

png

这个图里面倾斜的数据比之前看的要更清楚一些。

我们还可以观察拟合其他量度，最常用的还有均方误差（mean squared error，MSE），平均绝对误差（mean absolute deviation，MAD）。让我们用Python实现这两个量度。后面我们用scikit-learn内置的量度来评估回归模型的效果：

def MSE(target, predictions):
    squared_deviation = np.power(target - predictions, 2)
    return np.mean(squared_deviation)

MSE(boston.target, predictions)

21.897779217687496

def MAD(target, predictions):
    absolute_deviation = np.abs(target - predictions)
    return np.mean(absolute_deviation)

MAD(boston.target, predictions)

3.2729446379969396

How it works...

MSE的计算公式是：

$E(\hat y_t - y_i)^2$

计算预测值与实际值的差，平方之后再求平均值。这其实就是我们寻找最佳相关系数时是目标。高斯－马尔可夫定理（Gauss-Markov theorem）实际上已经证明了线性回归的回归系数的最佳线性无偏估计（BLUE）就是最小均方误差的无偏估计（条件是误差变量不相关，0均值，同方差）。在用岭回归弥补线性回归的不足主题中，我们会看到，当我们的相关系数是有偏估计时会发生什么。

MAD是平均绝对误差，计算公式为：

$E|\hat y_t - y_i|$

线性回归的时候MAD通常不用，但是值得一看。为什么呢？可以看到每个量度的情况，还可以判断哪个量度更重要。例如，用MSE，较大的误差会获得更大的惩罚，因为平方把它放大。

There's more...

还有一点需要说明，那就是相关系数是随机变量，因此它们是有分布的。让我们用bootstrapping（重复试验）来看看犯罪率的相关系数的分布情况。bootstrapping是一种学习参数估计不确定性的常用手段：

n_bootstraps = 1000
len_boston = len(boston.target)
subsample_size = np.int(0.5*len_boston)
subsample = lambda: np.random.choice(np.arange(0, len_boston),size=subsample_size)
coefs = np.ones(n_bootstraps) #相关系数初始值设为1
for i in range(n_bootstraps):
    subsample_idx = subsample()
    subsample_X = boston.data[subsample_idx]
    subsample_y = boston.target[subsample_idx]
    lr.fit(subsample_X, subsample_y)
    coefs[i] = lr.coef_[0]

我们可以看到这个相关系数的分布直方图：

f = plt.figure(figsize=(7, 5))
ax = f.add_subplot(111)
ax.hist(coefs, bins=50, color='b', alpha=.5)
ax.set_title("Histogram of the lr.coef_[0].");

png

我们还想看看重复试验后的置信区间：

np.percentile(coefs, [2.5, 97.5])

array([-0.18030624,  0.03816062])

置信区间的范围表面犯罪率其实不影响房价，因为0在置信区间里面，表面犯罪率可能与房价无关。

值得一提的是，bootstrapping可以获得更好的相关系数估计值，因为使用bootstrapping方法的均值，会比普通估计方法更快地**收敛（converge）**到真实均值。

2.3 用岭回归弥补线性回归的不足

本主题将介绍岭回归。和线性回归不同，它引入了正则化参数来“缩减”相关系数。当数据集中存在共线因素时，岭回归会很有用。

Getting ready

让我们加载一个不满秩（low effective rank）数据集来比较岭回归和线性回归。秩是矩阵线性无关组的数量，满秩是指一个 $m \times n$ 矩阵中行向量或列向量中现行无关组的数量等于 $min(m,n)$ 。

How to do it...

首先我们用make_regression建一个有3个自变量的数据集，但是其秩为2，因此3个自变量中有两个自变量存在相关性。

from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=2000, n_features=3, effective_rank=2, noise=10)

首先，我们用普通的线性回归拟合：

import numpy as np
from sklearn.linear_model import LinearRegression
lr = LinearRegression()

def fit_2_regression(lr):
    n_bootstraps = 1000
    coefs = np.ones((n_bootstraps, 3))
    len_data = len(reg_data)
    subsample_size = np.int(0.75*len_data)
    subsample = lambda: np.random.choice(np.arange(0, len_data), size=subsample_size)

    for i in range(n_bootstraps):
        subsample_idx = subsample()
        subsample_X = reg_data[subsample_idx]
        subsample_y = reg_target[subsample_idx]
        lr.fit(subsample_X, subsample_y)
        coefs[i][0] = lr.coef_[0]
        coefs[i][1] = lr.coef_[1]
        coefs[i][2] = lr.coef_[2]
    %matplotlib inline
    import matplotlib.pyplot as plt
    f, axes = plt.subplots(nrows=3, sharey=True, sharex=True, figsize=(7, 5))
    f.tight_layout()

    for i, ax in enumerate(axes):
        ax.hist(coefs[:, i], color='b', alpha=.5)
        ax.set_title("Coef {}".format(i))
    return coefs

coefs = fit_2_regression(lr)

png

我们再用Ridge来拟合数据，对比结果：

from sklearn.linear_model import Ridge
coefs_r = fit_2_regression(Ridge())

png

两个回归算法的结果看着好像差不多，其实不然。岭回归的系数更接近0。让我们看看两者系数的差异：

np.mean(coefs - coefs_r, axis=0)

array([ 30.54319761,  25.1726559 ,   7.40345307])

从均值上看，线性回归比岭回归的系数要大很多。均值显示的差异其实是线性回归的系数隐含的偏差。那么，岭回归究竟有什么好处呢？让我们再看看系数的方差：

np.var(coefs, axis=0)

array([ 302.16242654,  177.36842779,  179.33610289])

np.var(coefs_r, axis=0)

array([ 19.60727206,  25.4807605 ,  22.74202917])

岭回归的系数方差也会小很多。这就是机器学习里著名的偏差-方差均衡(Bias-Variance Trade-off)。下一个主题我们将介绍如何调整岭回归的参数正则化，那是偏差-方差均衡的核心内容。

How it works...

介绍参数正则化之前，我们总结一下岭回归与线性回归的不同。前面介绍过，线性回归的目标是最小化 ${\begin{Vmatrix} \hat y - X \beta \end{Vmatrix}}^2$ 。

岭回归的目标是最小化 ${\begin{Vmatrix} \hat y - X \beta \end{Vmatrix}}^2 + {\begin{Vmatrix} \Gamma X \end{Vmatrix}}^2$ 。

其中， $\Gamma$ 就是岭回归Ridge的alpha参数，指单位矩阵的倍数。上面的例子用的是默认值。我们可以看看岭回归参数：

Ridge()

Ridge(alpha=1.0, copy_X=True, fit_intercept=True, max_iter=None,
   normalize=False, solver='auto', tol=0.001)

岭回归系数的解是： $\beta = {(X^TX + \Gamma ^ T \Gamma)}^{-1}X \hat y$

前面的一半和线性回归的系数的解是一样的，多了 $\Gamma ^ T \Gamma)$ 一项。矩阵 $A$ 的 $AA^T$ 的结果是对称矩阵，且是半正定矩阵（对任意非0向量 $x$ ，有 $x^TAx \ge 0$ ）。相当于在线性回归的目标函数分母部分增加了一个很大的数。这样就把系数挤向0了。这样的解释比较粗糙，要深入了解，建议你看看SVD（矩阵奇异值分解）与岭回归的关系。

2.4 优化岭回归参数

当你使用岭回归模型进行建模时，需要考虑Ridge的alpha参数。

例如，用OLS（普通最小二乘法）做回归也许可以显示两个变量之间的某些关系；但是，当alpha参数正则化之后，那些关系就会消失。做决策时，这些关系是否需要考虑就显得很重要了。

Getting ready

这是我们第一个进行模型参数优化的主题，通常用交叉检验（cross validation）完成。在后面的主题中，还会有更简便的方式实现这些，但是这里我们一步一步来实现岭回归的优化。

在scikit-learn里面，岭回归的 $\Gamma$ 参数就是RidgeRegression的alpha参数；因此，问题就是最优的alpha参数是什么。首先我们建立回归数据集：

from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=100, n_features=2, effective_rank=1, noise=10)

How to do it...

在linear_models模块中，有一个对象叫RidgeCV，表示岭回归交叉检验（ridge cross-validation）。这个交叉检验类似于留一交叉验证法（leave-one-out cross-validation，LOOCV）。这种方法是指训练数据时留一个样本，测试的时候用这个未被训练过的样本：

import numpy as np
from sklearn.linear_model import RidgeCV
rcv = RidgeCV(alphas=np.array([.1, .2, .3, .4]))
rcv.fit(reg_data, reg_target)

RidgeCV(alphas=array([ 0.1,  0.2,  0.3,  0.4]), cv=None, fit_intercept=True,
    gcv_mode=None, normalize=False, scoring=None, store_cv_values=False)

拟合模型之后，alpha参数就是最优参数：

rcv.alpha_

0.10000000000000001

这里，0.1是最优参数，我们还想看到0.1附近更精确的值：

rcv = RidgeCV(alphas=np.array([.08, .09, .1, .11, .12]))
rcv.fit(reg_data, reg_target)

RidgeCV(alphas=array([ 0.08,  0.09,  0.1 ,  0.11,  0.12]), cv=None,
    fit_intercept=True, gcv_mode=None, normalize=False, scoring=None,
    store_cv_values=False)

rcv.alpha_

0.080000000000000002

可以按照这个思路一直优化下去，这里只做演示，后面还是介绍更好的方法。

How it works...

上面的演示很直接，但是我们介绍一下为什么这么做，以及哪个值才是最优的。在交叉检验的每一步里，模型的拟合效果都是用测试样本的误差表示。默认情况使用平方误差。更多细节见There's more...一节。

我们可以让RidgeCV储存交叉检验的数据，这样就可以可视化整个过程：

alphas_to_test = np.linspace(0.0001, 0.05)
rcv3 = RidgeCV(alphas=alphas_to_test, store_cv_values=True)
rcv3.fit(reg_data, reg_target)

RidgeCV(alphas=array([ 0.0001 ,  0.00112,  0.00214,  0.00316,  0.00417,  0.00519,
        0.00621,  0.00723,  0.00825,  0.00927,  0.01028,  0.0113 ,
        0.01232,  0.01334,  0.01436,  0.01538,  0.01639,  0.01741,
        0.01843,  0.01945,  0.02047,  0.02149,  0.0225 ,  0.02352,
        0.02454,  0.02556...4185,
        0.04287,  0.04389,  0.04491,  0.04593,  0.04694,  0.04796,
        0.04898,  0.05   ]),
    cv=None, fit_intercept=True, gcv_mode=None, normalize=False,
    scoring=None, store_cv_values=True)

你会看到，我们测试了0.0001到0.05区间中的50个点。由于我们把store_cv_values设置成true，我们可以看到每一个值对应的拟合效果：

rcv3.cv_values_.shape

(100, 50)

通过100个样本的回归数据集，我们获得了50个不同的alpha值。我们可以看到50个误差值，最小的均值误差对应最优的alpha值：

smallest_idx = rcv3.cv_values_.mean(axis=0).argmin()
alphas_to_test[smallest_idx]

0.014357142857142857

此时问题转化成了“RidgeCV认可我们的选择吗？”可以再用下面的命令获取alpha值：

rcv3.alpha_

0.014357142857142857

通过可视化图形可以更直观的显示出来。我们画出50个测试alpha值的图：

%matplotlib inline
import matplotlib.pyplot as plt
f, ax = plt.subplots(figsize=(7, 5))
ax.set_title(r"Various values of $\alpha$")

xy = (alphas_to_test[smallest_idx], rcv3.cv_values_.mean(axis=0)[smallest_idx])
xytext = (xy[0] + .01, xy[1] + .1)

ax.annotate(r'Chosen $\alpha$', xy=xy, xytext=xytext,
            arrowprops=dict(facecolor='black', shrink=0, width=0)
            )
ax.plot(alphas_to_test, rcv3.cv_values_.mean(axis=0));

png

There's more...

如果我们想用其他误差自定义评分函数，也是可以实现的。前面我们介绍过MAD误差，我们可以用它来评分。首先我们需要定义损失函数：

def MAD(target, prediction):
    absolute_deviation = np.abs(target - prediction)
    return absolute_deviation.mean()

定义损失函数之后，我们用sklearn量度中的make_scorer函数来处理。这样做可以标准化自定义的函数，让scikit-learn对象可以使用它。另外，由于这是一个损失函数不是一个评分函数，是越低越好，所以要用sklearn来把最小化问题转化成最大化问题：

import sklearn
MAD = sklearn.metrics.make_scorer(MAD, greater_is_better=False)
rcv4 = RidgeCV(alphas=alphas_to_test, store_cv_values=True, scoring=MAD)
rcv4.fit(reg_data, reg_target)
smallest_idx = rcv4.cv_values_.mean(axis=0).argmin()
alphas_to_test[smallest_idx]

0.050000000000000003

2.5 LASSO正则化

LASSO（ least absolute shrinkage and selection operator，最小绝对值收缩和选择算子）方法与岭回归和LARS（least angle regression，最小角回归）很类似。与岭回归类似，它也是通过增加惩罚函数来判断、消除特征间的共线性。与LARS相似的是它也可以用作参数选择，通常得出一个相关系数的稀疏向量。

Getting ready

岭回归也不是万能药。有时就需要用LASSO回归来建模。本主题将用不同的损失函数，因此就要用对应的效果评估方法。

How to do it...

首先，我们还是用make_regression函数来建立数据集：

from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=200, n_features=500, n_informative=5, noise=5)

之后，我们导入lasso对象：

from sklearn.linear_model import Lasso
lasso = Lasso()

lasso包含很多参数，但是最意思的参数是alpha，用来调整 lasso的惩罚项，在**How it works...**会具体介绍。现在我们用默认值1。另外，和岭回归类似，如果设置为0，那么lasso就是线性回归：

lasso.fit(reg_data, reg_target)

Lasso(alpha=1.0, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)

再让我们看看还有多少相关系数非零：

import numpy as np
np.sum(lasso.coef_ != 0)

lasso_0 = Lasso(0)
lasso_0.fit(reg_data, reg_target)
np.sum(lasso_0.coef_ != 0)

d:\programfiles\Miniconda3\lib\site-packages\IPython\kernel\__main__.py:2: UserWarning: With alpha=0, this algorithm does not converge well. You are advised to use the LinearRegression estimator
  from IPython.kernel.zmq import kernelapp as app
d:\programfiles\Miniconda3\lib\site-packages\sklearn\linear_model\coordinate_descent.py:432: UserWarning: Coordinate descent with alpha=0 may lead to unexpected results and is discouraged.
  positive)





500

和我们设想的一样，如果用线性回归，没有一个相关系数变成0。而且，如果你这么运行代码，scikit-learn会给出建议，就像上面显示的那样。

How it works...

对线性回归来说，我们是最小化残差平方和。而LASSO回归里，我们还是最小化残差平方和，但是加了一个惩罚项会导致稀疏。如下所示：

{\begin{Vmatrix} \beta \end{Vmatrix}}_1}

最小化残差平方和的另一种表达方式是：

{\begin{Vmatrix} \beta \end{Vmatrix}}_1 \lt \beta

这个约束会让数据稀疏。LASSO回归的约束创建了围绕原点的超立方体（相关系数是轴），也就意味着大多数点都在各个顶点上，那里相关系数为0。而岭回归创建的是超平面，因为其约束是L2范数，少一个约束，但是即使有限制相关系数也不会变成0。

LASSO交叉检验

上面的公式中，选择适当的 $\lambda$ （在scikit-learn的Lasso里面是alpha，但是书上都是 $\lambda$ ）参数是关键。我们可以自己设置，也可以通过交叉检验来获取最优参数：

from sklearn.linear_model import LassoCV
lassocv = LassoCV()
lassocv.fit(reg_data, reg_target)

LassoCV(alphas=None, copy_X=True, cv=None, eps=0.001, fit_intercept=True,
    max_iter=1000, n_alphas=100, n_jobs=1, normalize=False, positive=False,
    precompute='auto', random_state=None, selection='cyclic', tol=0.0001,
    verbose=False)

lassocv有一个属性就是确定最合适的 $\lambda$ ：

lassocv.alpha_

0.58535963603062136

计算的相关系数也可以看到：

lassocv.coef_[:5]

array([ 0.       , -0.       ,  0.       ,  0.0192606, -0.       ])

用最近的参数拟合后，lassocv的非零相关系数有29个：

np.sum(lassocv.coef_ != 0)

LASSO特征选择

LASSO通常用来为其他方法所特征选择。例如，你可能会用LASSO回归获取适当的特征变量，然后在其他算法中使用。

要获取想要的特征，需要创建一个非零相关系数的列向量，然后再其他算法拟合：

mask = lassocv.coef_ != 0
new_reg_data = reg_data[:, mask]
new_reg_data.shape

(200, 29)

2.6 LARS正则化

如果斯坦福大学的Bradley Efron, Trevor Hastie, Iain Johnstone和Robert Tibshirani没有发现它的话[1]，LARS(Least Angle Regression，最小角回归)可能有一天会被你想出来，它借用了威廉·吉尔伯特·斯特朗（William Gilbert Strang）介绍过的高斯消元法（Gaussian elimination）的灵感。

Getting ready

LARS是一种回归手段，适用于解决高维问题，也就是 $p >> n$ 的情况，其中 $p$ 表示列或者特征变量， $n$ 表示样本数量。

How to do it...

首先让我们导入必要的对象。这里我们用的数据集是200个数据，500个特征。我们还设置了一个低噪声，和少量提供信息的（informative）特征：

import numpy as np
from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=200,n_features=500, n_informative=10, noise=2)

由于我们用了10个信息特征，因此我们还要为LARS设置10个非0的相关系数。我们事先可能不知道信息特征的准确数量，但是出于试验的目的是可行的：

from sklearn.linear_model import Lars
lars = Lars(n_nonzero_coefs=10)
lars.fit(reg_data, reg_target)

Lars(copy_X=True, eps=2.2204460492503131e-16, fit_intercept=True,
   fit_path=True, n_nonzero_coefs=10, normalize=True, precompute='auto',
   verbose=False)

我们可以检验一下看看LARS的非0相关系数的和：

np.sum(lars.coef_ != 0)

问题在于为什么少量的特征反而变得更加有效。要证明这一点，让我们用一半数量来训练两个LARS模型，一个用12个非零相关系数，另一个非零相关系数用默认值。这里用12个是因为我们对重要特征的数量有个估计，但是可能无法确定准确的数量：

train_n = 100
lars_12 = Lars(n_nonzero_coefs=12)
lars_12.fit(reg_data[:train_n], reg_target[:train_n])

Lars(copy_X=True, eps=2.2204460492503131e-16, fit_intercept=True,
   fit_path=True, n_nonzero_coefs=12, normalize=True, precompute='auto',
   verbose=False)

lars_500 = Lars() #默认就是500
lars_500.fit(reg_data[:train_n], reg_target[:train_n])

Lars(copy_X=True, eps=2.2204460492503131e-16, fit_intercept=True,
   fit_path=True, n_nonzero_coefs=500, normalize=True, precompute='auto',
   verbose=False)

现在，让我们看看拟合数据的效果如何，如下所示：

np.mean(np.power(reg_target[train_n:] - lars.predict(reg_data[train_n:]), 2))

18.607806437043894

np.mean(np.power(reg_target[train_n:] - lars_12.predict(reg_data[train_n:]), 2))

529.97993250189643

np.mean(np.power(reg_target[train_n:] - lars_500.predict(reg_data[train_n:]), 2))

2.3236770314162846e+34

仔细看看这组结果；测试集的误差明显高很多。高维数据集问题就在于此；通常面对大量的特征时，想找出一个对训练集拟合很好的模型并不难，但是拟合过度却是更大的问题。

How it works...

LARS通过重复选择与残存变化相关的特征。从图上看，相关性实际上就是特征与残差之间的最小角度；这就是LARS名称的由来。

选择第一个特征之后，LARS会继续沿着最小角的方向移动，直到另一个特征与残差有同样数量的相关性。然后，LARS会沿着两个特征组合的角度移动。如下图所示：

%matplotlib inline
import matplotlib.pyplot as plt
def unit(*args):
    squared = map(lambda x: x**2, args)
    distance = sum(squared) ** (.5)
    return map(lambda x: x / distance, args)

f, ax = plt.subplots(nrows=3, figsize=(5, 10))
plt.tight_layout()
ax[0].set_ylim(0, 1.1)
ax[0].set_xlim(0, 1.1)

x, y = unit(1, 0.02)
ax[0].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[0].text(x + .05, y + .05, r"$x_1$")

x, y = unit(.5, 1)
ax[0].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[0].text(x + .05, y + .05, r"$x_2$")

x, y = unit(1, .45)
ax[0].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[0].text(x + .05, y + .05, r"$y$")

ax[0].set_title("No steps")

# step 1
ax[1].set_title("Step 1")
ax[1].set_ylim(0, 1.1)
ax[1].set_xlim(0, 1.1)

x, y = unit(1, 0.02)
ax[1].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[1].text(x + .05, y + .05, r"$x_1$")

x, y = unit(.5, 1)
ax[1].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[1].text(x + .05, y + .05, r"$x_2$")

x, y = unit(.5, 1)
ax[1].arrow(.5, 0.01, x, y, ls='dashed', edgecolor='black', facecolor='black')
ax[1].text(x + .5 + .05, y + .01 + .05, r"$x_2$")

ax[1].arrow(0, 0, .47, .01, width=.0015, edgecolor='black', facecolor='black')
ax[1].text(.47-.15, .01 + .03, "Step 1")

x, y = unit(1, .45)
ax[1].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[1].text(x + .05, y + .05, r"$y$")

# step 2
ax[2].set_title("Step 2")
ax[2].set_ylim(0, 1.1)
ax[2].set_xlim(0, 1.1)

x, y = unit(1, 0.02)
ax[2].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[2].text(x + .05, y + .05, r"$x_1$")

x, y = unit(.5, 1)
ax[2].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[2].text(x + .05, y + .05, r"$x_2$")

x, y = unit(.5, 1)
ax[2].arrow(.5, 0.01, x, y, ls='dashed', edgecolor='black', facecolor='black')
ax[2].text(x + .5 + .05, y + .01 + .05, r"$x_2$")

ax[2].arrow(0, 0, .47, .01, width=.0015, edgecolor='black', facecolor='black')
ax[2].text(.47-.15, .01 + .03, "Step 1")

##  step 2
x, y = unit(1, .45)
ax[2].arrow(.5, .02, .4, .35, width=.0015, edgecolor='black', facecolor='black')
ax[2].text(x, y - .1, "Step 2")

x, y = unit(1, .45)
ax[2].arrow(0, 0, x, y, edgecolor='black', facecolor='black')
ax[2].text(x + .05, y + .05, r"$y$");

png

具体过程是，我们把 $x2$ 沿着 $x1$ 方向移动到一个位置： $x1$ 与 $y$ 的点积与 $x1$ 与 $y$ 的点积相同。到了这个位置之后，我们再沿着 $x1$ 和 $x2$ 夹角的一半的方向移动。

There's more...

和我们前面用交叉检验来优化领回归模型一样，我们可以对LARS做交叉检验：

from sklearn.linear_model import LarsCV
lcv = LarsCV()
lcv.fit(reg_data, reg_target)

d:\Miniconda3\lib\site-packages\sklearn\linear_model\least_angle.py:285: ConvergenceWarning: Regressors in active set degenerate. Dropping a regressor, after 168 iterations, i.e. alpha=2.278e-02, with an active set of 132 regressors, and the smallest cholesky pivot element being 6.144e-08
  ConvergenceWarning)
d:\Miniconda3\lib\site-packages\sklearn\linear_model\least_angle.py:285: ConvergenceWarning: Regressors in active set degenerate. Dropping a regressor, after 168 iterations, i.e. alpha=2.105e-02, with an active set of 132 regressors, and the smallest cholesky pivot element being 9.771e-08
  ConvergenceWarning)





LarsCV(copy_X=True, cv=None, eps=2.2204460492503131e-16, fit_intercept=True,
    max_iter=500, max_n_alphas=1000, n_jobs=1, normalize=True,
    precompute='auto', verbose=False)

用交叉检验可以帮助我们确定需要使用的非零相关系数的最佳数量。验证如下所示：

np.sum(lcv.coef_ != 0)

说实话，LARS的精髓还没有领会，抽空会把原文译出来，看各种解释不如看原文。

[1] Efron, Bradley; Hastie, Trevor; Johnstone, Iain and Tibshirani, Robert(2004). "Least Angle Regression". Annals of Statistics 32(2): pp. 407–499.doi:10.1214/009053604000000067. MR 2060166.

2.7 用线性方法处理分类问题——逻辑回归

实际上线性模型也可以用于分类任务。方法是把一个线性模型拟合成某个类型的概率分布，然后用一个函数建立阈值来确定结果属于哪一类。

Getting ready

这里用的函数是经典的逻辑函数。一个非常简单的函数：

$f(x)= \frac 1 {1+e^{-t}}$

它的图形如下图所示：

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
f, ax = plt.subplots(figsize=(10, 5))
rng = np.linspace(-5, 5)
log_f = np.apply_along_axis(lambda x:1 / (1 + np.exp(-x)), 0, rng)
ax.set_title("Logistic Function between [-5, 5]")
ax.plot(rng, log_f);

png

让我们用make_classification方法创建一个数据集来进行分类：

from sklearn.datasets import make_classification
X, y = make_classification(n_samples=1000, n_features=4)

How to do it...

LogisticRegression对象和其他线性模型的用法一样：

from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()

我们将把前面200个数据作为训练集，最后200个数据作为测试集。因为这是随机数据集，所以用最后200个数据没问题。但是如果处理具有某种结构的数据，就不能这么做了（例如，你的数据集是时间序列数据）：

X_train = X[:-200]
X_test = X[-200:]
y_train = y[:-200]
y_test = y[-200:]

在本书后面的内容里，我们将详细介绍交叉检验。这里，我们需要的只是用逻辑回归拟合模型。我们会关注训练集的预测结果，就像测试集预测结果一样。经常对比两个数据集预测正确率是个好做法。通常，你在训练集获得的结果更好；模型在测试集上预测失败的比例也至关重要：

lr.fit(X_train, y_train)
y_train_predictions = lr.predict(X_train)
y_test_predictions = lr.predict(X_test)

现在我们有了预测值，让我们看看预测的效果。这里，我们只简单看看预测正确的比例；后面，我们会详细的介绍分类模型效果的评估方法。

计算很简单，就是用预测正确的数量除以总样本数：

(y_train_predictions == y_train).sum().astype(float) / y_train.shape[0]

0.89375000000000004

测试集的效果是：

(y_test_predictions == y_test).sum().astype(float) / y_test.shape[0]

0.90500000000000003

可以看到，测试集的正确率和训练集的结果差不多。但是实际中通常差别很大。

现在问题变成，怎么把逻辑函数转换成分类方法。

首先，线性回归希望找到一个线性方程拟合出给定自变量 $X$ 条件下因变量 $Y$ 的期望值，就是 $E(Y|X)=x \beta$ 。这里 $Y$ 的值是某个类型发生的概率。因此，我们要解决的分类问题就是 $E(p|X)=x \beta$ 。然后，只要阈值确定，就会有 $Logit(p) = X \beta$ 。这个理念的扩展形式可以构成许多形式的回归行为，例如，泊松过程（Poisson）。

There's more...

下面的内容你以后肯定会遇到。一种情况是一个类型与其他类型的权重不同；例如，一个能可能权重很大，99%。这种情况在分类工作中经常遇到。经典案例就是信用卡虚假交易检测，大多数交易都不是虚假交易，但是不同类型误判的成本相差很大。

让我们建立一个分类问题，类型 $y$ 的不平衡权重95%，我们看看基本的逻辑回归模型如何处理这类问题：

X, y = make_classification(n_samples=5000, n_features=4, weights=[.95])
sum(y) / (len(y)*1.) #检查不平衡的类型

0.0562

建立训练集和测试集，然后用逻辑回归拟合：

X_train = X[:-500]
X_test = X[-500:]
y_train = y[:-500]
y_test = y[-500:]

lr.fit(X_train, y_train)
y_train_predictions = lr.predict(X_train)
y_test_predictions = lr.predict(X_test)

现在我们在看看模型拟合的情况：

(y_train_predictions == y_train).sum().astype(float) / y_train.shape[0]

0.96711111111111114

(y_test_predictions == y_test).sum().astype(float) / y_test.shape[0]

0.96799999999999997

结果看着还不错，但这是说如果我们把一个交易预测成正常交易（或者称为类型0），那么我们有95%左右的可能猜对。如果我们想看看模型对类型1的预测情况，可能就不是那么好了：

(y_test[y_test==1] == y_test_predictions[y_test==1]).sum().astype(float) / y_test[y_test==1].shape[0]

0.5

如果相比正常交易，我们更关心虚假交易；那么这是由商业规则决定的，我们可能会改变预测正确和预测错误的权重。

通常情况下，虚假交易与正常交易的权重与训练集的类型权重的倒数一致。但是，因为我们更关心虚假交易，所有让我们用多重采样（oversample）方法来表示虚假交易与正常交易的权重。

lr = LogisticRegression(class_weight={0: .15, 1: .85})
lr.fit(X_train, y_train)

LogisticRegression(C=1.0, class_weight={0: 0.15, 1: 0.85}, dual=False,
          fit_intercept=True, intercept_scaling=1, max_iter=100,
          multi_class='ovr', penalty='l2', random_state=None,
          solver='liblinear', tol=0.0001, verbose=0)

让我们再预测一下结果：

y_train_predictions = lr.predict(X_train)
y_test_predictions = lr.predict(X_test)

(y_test[y_test==1] == y_test_predictions[y_test==1]).sum().astype(float) / y_test[y_test==1].shape[0]

0.70833333333333337

但是，这么做需要付出什么代价？让我们看看：

(y_test_predictions == y_test).sum().astype(float) / y_test.shape[0]

0.93999999999999995

可以看到，准确率降低了3%。这样是否合适由你的问题决定。如果与虚假交易相关的评估成本非常高，那么它就能抵消追踪虚假交易付出的成本。

2.8 贝叶斯岭回归

在用岭回归弥补线性回归的不足主题中，我们介绍了岭回归优化的限制条件。我们还介绍了相关系数的先验概率分布的贝叶斯解释，将很大程度地影响着先验概率分布，先验概率分布通常均值是0。

因此，现在我们就来演示如何scikit-learn来应用这种解释。

Getting ready

岭回归和套索回归（lasso regression）用贝叶斯观点来解释，与频率优化观点解释相反。scikit-learn只实现了贝叶斯岭回归，但是在*How it works...*一节，我们将对比两种回归算法。

首先，我们创建一个回归数据集：

from sklearn.datasets import make_regression
X, y = make_regression(1000, 10, n_informative=2, noise=20)

How to do it...

我们可以把岭回归加载进来拟合模型：

from sklearn.linear_model import BayesianRidge
br = BayesianRidge()

有两组相关系数，分别是alpha_1 / alpha_2和lambda_1 / lambda_2。其中，alpha_*是先验概率分布的 $\alpha$ 超参数，lambda_*是先验概率分布的 $\lambda$ 超参数。

首先，让我们不调整参数直接拟合模型：

br.fit(X, y)
br.coef_

array([ -1.39241213,   0.14671513,  -0.08150797,  37.50250891,
         0.21850082,  -0.78482779,  -0.26717555,  -0.71319956,
         0.7926308 ,   5.74658302])

现在，我们来调整超参数，注意观察相关系数的变化：

br_alphas = BayesianRidge(alpha_1=10, lambda_1=10)
br_alphas.fit(X, y)
br_alphas.coef_

array([ -1.38807423,   0.14050794,  -0.08309391,  37.3032803 ,
         0.2254332 ,  -0.77031801,  -0.27005478,  -0.71632657,
         0.78501276,   5.71928608])

How it works...

因为是贝叶斯岭回归，我们假设先验概率分布带有误差和 $\alpha$ 参数，先验概率分布都服从 $\Gamma$ 分布。

$\Gamma$ 分布是一种极具灵活性的分布。不同的形状参数和尺度参数的 $\Gamma$ 分布形状有差异。1e-06是 scikit-learn里面BayesianRidge形状参数的默认参数值。

import matplotlib.pyplot as plt
%matplotlib inline

from scipy.stats import gamma
import numpy as np

form = lambda x, y: "loc={}, scale={}".format(x, y)
g = lambda x, y=1e-06, z=1e-06: gamma.pdf(x, y, z)
g2 = lambda x, y=1e-06, z=1: gamma.pdf(x, y, z)
g3 = lambda x, y=1e-06, z=2: gamma.pdf(x, y, z)
rng = np.linspace(0, 5)
f, ax = plt.subplots(figsize=(8, 5))

ax.plot(rng, list(map(g, rng)), label=form(1e-06, 1e-06), color='r')
ax.plot(rng, list(map(g2, rng)), label=form(1e-06, 1), color='g')
ax.plot(rng, list(map(g3, rng)), label=form(1e-06, 2), color='b')
ax.set_title("Different Shapes of the Gamma Distribution")
ax.legend();

png

你会看到，相关系数最终都会收缩到0，尤其当形状参数特别小的时候。

There's more...

就像我前面介绍的，还有一种套索回归的贝叶斯解释。我们把先验概率分布看出是相关系数的函数；它们本身都是随机数。对于套索回归，我们选择一个可以产生0的分布，比如双指数分布（Double Exponential Distribution，也叫Laplace distribution）。

from scipy.stats import laplace
form = lambda x, y: "loc={}, scale={}".format(x, y)
g = lambda x: laplace.pdf(x)
rng = np.linspace(-5, 5)
f, ax = plt.subplots(figsize=(8, 5))

ax.plot(rng, list(map(g, rng)), color='r')
ax.set_title("Example of Double Exponential Distribution");

png

留意看x轴为0处的顶点。这将会使套索回归的相关系数为0。通过调整超参数，还有可能创建出相关系数为0的情况，这由问题的具体情况决定。

2.9 用梯度提升回归从误差中学习

梯度提升回归（Gradient boosting regression，GBR）是一种从它的错误中进行学习的技术。它本质上就是集思广益，集成一堆较差的学习算法进行学习。有两点需要注意：

每个学习算法准备率都不高，但是它们集成起来可以获得很好的准确率。
这些学习算法依次应用，也就是说每个学习算法都是在前一个学习算法的错误中学习

Getting ready

我们还是用基本的回归数据来演示GBR：

import numpy as np
from sklearn.datasets import make_regression
X, y = make_regression(1000, 2, noise=10)

How to do it...

GBR算是一种集成模型因为它是一个集成学习算法。这种称谓的含义是指GBR用许多较差的学习算法组成了一个更强大的学习算法：

from sklearn.ensemble import GradientBoostingRegressor as GBR
gbr = GBR()
gbr.fit(X, y)
gbr_preds = gbr.predict(X)

很明显，这里应该不止一个模型，但是这种模式现在很简明。现在，让我们用基本回归算法来拟合数据当作参照：

from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X, y)
lr_preds = lr.predict(X)

有了参照之后，让我们看看GBR算法与线性回归算法效果的对比情况。图像生成可以参照第一章正态随机过程的相关主题，首先需要下面的计算：

gbr_residuals = y - gbr_preds
lr_residuals = y - lr_preds

%matplotlib inline
from matplotlib import pyplot as plt

f, ax = plt.subplots(figsize=(7, 5))
f.tight_layout()
ax.hist(gbr_residuals,bins=20,label='GBR Residuals', color='b', alpha=.5);
ax.hist(lr_residuals,bins=20,label='LR Residuals', color='r', alpha=.5);
ax.set_title("GBR Residuals vs LR Residuals")
ax.legend(loc='best');

png

看起来好像GBR拟合的更好，但是并不明显。让我们用95%置信区间（Confidence interval,CI）对比一下：

np.percentile(gbr_residuals, [2.5, 97.5])

array([-16.73937398,  15.96258406])

np.percentile(lr_residuals, [2.5, 97.5])

array([-19.03378242,  19.45950191])

GBR的置信区间更小，数据更集中，因此其拟合效果更好；我们还可以对GBR算法进行一些调整来改善效果。我用下面的例子演示一下，然后在下一节介绍优化方法：

n_estimators = np.arange(100, 1100, 350)
gbrs = [GBR(n_estimators=n_estimator) for n_estimator in n_estimators]
residuals = {}
for i, gbr in enumerate(gbrs):
    gbr.fit(X, y)
    residuals[gbr.n_estimators] = y - gbr.predict(X)

f, ax = plt.subplots(figsize=(7, 5))
f.tight_layout()
colors = {800:'r', 450:'g', 100:'b'}
for k, v in residuals.items():
    ax.hist(v,bins=20,label='n_estimators: %d' % k, color=colors[k], alpha=.5);
ax.set_title("Residuals at Various Numbers of Estimators")
ax.legend(loc='best');

png

图像看着有点混乱，但是依然可以看出随着估计器数据的增加，误差在减少。不过，这并不是一成不变的。首先，我们没有交叉检验过，其次，随着估计器数量的增加，训练时间也会变长。现在我们用数据比较小没什么关系，但是如果数据再放大一两倍问题就出来了。

How it works...

上面例子中GBR的第一个参数是n_estimators，指GBR使用的学习算法的数量。通常，如果你的设备性能更好，可以把n_estimators设置的更大，效果也会更好。还有另外几个参数要说明一下。

你应该在优化其他参数之前先调整max_depth参数。因为每个学习算法都是一颗决策树，max_depth决定了树生成的节点数。选择合适的节点数量可以更好的拟合数据，而更多的节点数可能造成拟合过度。

loss参数决定损失函数，也直接影响误差。ls是默认值，表示最小二乘法（least squares）。还有最小绝对值差值，Huber损失和分位数损失（quantiles）等等。

第三章使用距离向量构建模型

作者：Trent Hauck

译者：飞龙

协议：CC BY-NC-SA 4.0

这一章中，我们会涉及到聚类。聚类通常和非监督技巧组合到一起。这些技巧假设我们不知道结果变量。这会使结果模糊，以及实践客观。但是，聚类十分有用。我们会看到，我们可以使用聚类，将我们的估计在监督设置中“本地化”。这可能就是聚类非常高效的原因。它可以处理很大范围的情况，通常，结果也不怎么正常。

这一章中我们会浏览大量应用，从图像处理到回归以及离群点检测。通过这些应用，我们会看到聚类通常可以通过概率或者优化结构来观察。不同的解释会导致不同的权衡。我们会看到，如何训练模型，以便让工具尝试不同模型，在面对聚类问题的时候。

3.1 使用 KMeans 对数据聚类

聚类是个非常实用的技巧。通常，我们在采取行动时需要分治。考虑公司的潜在客户列表。公司可能需要将客户按类型分组，之后为这些分组划分职责。聚类可以使这个过程变得容易。

KMeans 可能是最知名的聚类算法之一，并且也是最知名的无监督学习技巧之一。

准备

首先，让我们看一个非常简单的聚类，之后我们再讨论 KMeans 如何工作。

>>> from sklearn.datasets import make_blobs 
>>> blobs, classes = make_blobs(500, centers=3)

同样，由于我们绘制一些图表，导入matplotlib，像这样：

>>> import matplotlib.pyplot as plt

操作步骤

我们打算浏览一个简单的例子，它对伪造数据进行聚类。之后我们会稍微谈论一下，KMeans 如何工作，来寻找最优的块数量。

看一看我们的数据块，我们可以看到，有三个不同的簇。

>>> f, ax = plt.subplots(figsize=(7.5, 7.5)) 
>>> ax.scatter(blobs[:, 0], blobs[:, 1], color=rgb[classes]) 
>>> rgb = np.array(['r', 'g', 'b']) 
>>> ax.set_title("Blobs")

输出如下：

现在我们可以使用 KMeans 来寻找这些簇的形心。第一个例子中，我们假装知道有三个形心。

>>> from sklearn.cluster import KMeans 
>>> kmean = KMeans(n_clusters=3) 
>>> kmean.fit(blobs) 
KMeans(copy_x=True, init='k-means++', max_iter=300, n_clusters=3,
       n_init=10, n_jobs=1, precompute_distances=True,
       random_state=None, tol=0.0001, verbose=0)

>>> kmean.cluster_centers_ array([[ 0.47819567,  1.80819197],

[ 0.08627847,  8.24102715], 
[ 5.2026125 ,  7.86881767]])

>>> f, ax = plt.subplots(figsize=(7.5, 7.5)) 
>>> ax.scatter(blobs[:, 0], blobs[:, 1], color=rgb[classes]) 
>>> ax.scatter(kmean.cluster_centers_[:, 0],
               kmean.cluster_centers_[:, 1], marker='*', s=250,
               color='black', label='Centers')
>>> ax.set_title("Blobs") 
>>> ax.legend(loc='best')

下面的截图展示了输出：

其它属性也很实用。例如，labels_属性会产生每个点的预期标签。

>>> kmean.labels_[:5] 
array([1, 1, 2, 2, 1], dtype=int32)

我们可以检查，例如，labels_是否和类别相同，但是由于 KMeans 不知道类别是什么，它不能给两个类别分配相同的索引值：

>>> classes[:5] 
array([0, 0, 2, 2, 0])

将类别中的1变成0来查看是否与labels_匹配。

transform函数十分有用，它会输出每个点到形心的距离。

>>> kmean.transform(blobs)[:5] 
array([[ 6.47297373,  1.39043536,  6.4936008 ],
       [ 6.78947843,  1.51914705,  3.67659072],
       [ 7.24414567,  5.42840092,  0.76940367],
       [ 8.56306214,  5.78156881,  0.89062961],
       [ 7.32149254,  0.89737788,  5.12246797]])

工作原理

KMeans 实际上是个非常简单的算法，它使簇中的点到均值的距离的平方和最小。

首先它会设置一个预定义的簇数量K，之后执行这些事情：

将每个数据点分配到最近的簇中。
通过计算初中每个数据点的均值，更新每个形心。

直到满足特定条件。

3.2 优化形心数量

形心难以解释，并且也难以判断是否数量正确。理解你的数据是否是未分类的十分重要，因为这会直接影响我们可用的评估手段。

准备

为无监督学习评估模型表现是个挑战。所以，在了解真实情况的时候，sklearn拥有多种方式来评估聚类，但在不了解时就很少。

我们会以一个简单的簇模型开始，并评估它的相似性。这更多是出于机制的目的，因为测量一个簇的相似性在寻找簇数量的真实情况时显然没有用。

操作步骤

为了开始，我们会创建多个数据块，它们可用于模拟数据簇。

>>> from sklearn.datasets import make_blobs 
>>> import numpy as np 
>>> blobs, classes = make_blobs(500, centers=3)

>>> from sklearn.cluster import KMeans 
>>> kmean = KMeans(n_clusters=3) 
>>> kmean.fit(blobs) 
KMeans(copy_x=True, init='k-means++', max_iter=300, n_clusters=3,
       n_init=10, n_jobs=1, precompute_distances=True,
       random_state=None, tol=0.0001, verbose=0)

首先，我们查看轮廓（Silhouette）距离。轮廓距离是簇内不相似性、最近的簇间不相似性、以及这两个值最大值的比值。它可以看做簇间分离程度的度量。

让我们看一看数据点到形心的距离分布，理解轮廓距离非常有用。

>>> from sklearn import metrics 
>>> silhouette_samples = metrics.silhouette_samples(blobs,
                         kmean.labels_) 
>>> np.column_stack((classes[:5], silhouette_samples[:5]))

array([[ 1.,  0.87617292],
       [ 1.,  0.89082363],
       [ 1.,  0.88544994],
       [ 1.,  0.91478369],
       [ 1.,  0.91308287]]) 
>>> f, ax = plt.subplots(figsize=(10, 5))

>>> ax.set_title("Hist of Silhouette Samples") 
>>> ax.hist(silhouette_samples)

输出如下：

要注意，通常接近 1 的系数越高，分数就越高。

工作原理

轮廓系数的均值通常用于描述整个模型的拟合度。

>>> silhouette_samples.mean() 
0.57130462953339578

这十分普遍，事实上，metrics模块提供了一个函数来获得刚才的值。

现在，让我们训练多个簇的模型，并看看平均得分是什么样：

# first new ground truth 
>>> blobs, classes = make_blobs(500, centers=10) 
>>> sillhouette_avgs = []
# this could take a while 
>>> for k in range(2, 60):
       kmean = KMeans(n_clusters=k).fit(blobs)
       sillhouette_avgs.append(metrics.silhouette_score(blobs,
                               kmean.labels_))
>>> f, ax = plt.subplots(figsize=(7, 5)) 
>>> ax.plot(sillhouette_avgs)

下面是输出：

这个绘图表明，轮廓均值随着形心数量的变化情况。我们可以看到最优的数量是 3，根据所生成的数据。但是最优的数量看起来是 6 或者 7。这就是聚类的实际情况，十分普遍，我们不能获得正确的簇数量，我们只能估计簇数量的近似值。

3.3 评估聚类的正确性

我们之前讨论了不知道真实情况的条件下的聚类评估。但是，我们还没有讨论簇已知条件下的 KMeans 评估。在许多情况下，这都是不可知的，但是如果存在外部的标注，我们就会知道真实情况，或者至少是代理。

准备

所以，让我们假设有一个世界，其中我们有一些外部代理，向我们提供了真实情况。

我们会创建一个简单的数据集，使用多种方式评估相对于真实庆康的正确性。之后讨论它们。

操作步骤

在我们开始度量之前，让我们先查看数据集：

>>> f, ax = plt.subplots(figsize=(7, 5))

>>> colors = ['r', 'g', 'b']

>>> for i in range(3):
       p = blobs[ground_truth == i]
       ax.scatter(p[:,0], p[:,1], c=colors[i],
       label="Cluster {}".format(i))
       
>>> ax.set_title("Cluster With Ground Truth") 
>>> ax.legend()

>>> f.savefig("9485OS_03-16")

下面是输出：

既然我们已经训练了模型，让我们看看簇的形心：

>>> f, ax = plt.subplots(figsize=(7, 5))

>>> colors = ['r', 'g', 'b']

>>> for i in range(3):
       p = blobs[ground_truth == i]
       ax.scatter(p[:,0], p[:,1], c=colors[i],        label="Cluster {}".format(i))
>>> ax.scatter(kmeans.cluster_centers_[:, 0],
               kmeans.cluster_centers_[:, 1], s=100,
               color='black',
               label='Centers') 
>>> ax.set_title("Cluster With Ground Truth") 
>>> ax.legend()
>>> f.savefig("9485OS_03-17")

下面是输出：

既然我们能够将聚类表现看做分类练习，在其语境中有用的方法在这里也有用：

>>> for i in range(3):
       print (kmeans.labels_ == ground_truth)[ground_truth == i]
       .astype(int).mean()
       
0.0778443113772 
0.990990990991 
0.0570570570571

很显然我们有一些错乱的簇。所以让我们将其捋直，之后我们查看准确度。

>>> new_ground_truth = ground_truth.copy()
>>> new_ground_truth[ground_truth == 0] = 2 
>>> new_ground_truth[ground_truth == 2] = 0

>>> for i in range(3):
       print (kmeans.labels_ == new_ground_truth)[ground_truth == i]
       .astype(int).mean()
       
0.919161676647 
0.990990990991 
0.90990990991

所以我们 90% 的情况下都是正确的。第二个相似性度量是互信息（ mutual information score）得分。

>>> from sklearn import metrics

>>> metrics.normalized_mutual_info_score(ground_truth, kmeans.labels_)

0.78533737204433651

分数靠近 0，就说明标签的分配可能不是按照相似过程生成的。但是分数靠近 1，就说明两个标签有很强的一致性。

例如，让我们看一看互信息分数自身的情况：

>>> metrics.normalized_mutual_info_score(ground_truth, ground_truth)

1.0

通过名称，我们可以分辨出可能存在未规范化的 mutual_info_score：

>>> metrics.mutual_info_score(ground_truth, kmeans.labels_)

0.78945287371677486

这非常接近了。但是，规范化的互信息是互信息除以每个真实值和标签的熵的乘积的平方根。

有一个度量方式我们尚未讨论，并且不依赖于真实情况，就是惯性（inertia）度量。当前，它作为一种度量并没有详细记录。但是，它是 KMeans 中最简单的度量。

惯性是每个数据点和它所分配的簇的平方差之和。我们可以稍微使用 NumPy 来计算它：

>>> kmeans.inertia_

3.4 使用 MiniBatch KMeans 处理更多数据

KMeans 是一个不错的方法，但是不适用于大量数据。这是因为 KMenas 的复杂度。也就是说，我们可以使用更低的算法复杂度来获得近似解。

准备

MiniBatch Kmeans 是 KMeans 的更快实现。KMeans 的计算量非常大，问题是 NPH 的。

但是，使用 MiniBatch KMeans，我们可以将 KMeans 加速几个数量级。这通过处理多个子样本来完成，它们叫做 MiniBatch。如果子样本是收敛的，并且拥有良好的初始条件，就得到了常规 KMeans 的近似解。

操作步骤

让我们对 MiniBatch 聚类做一个概要的性能分析。首先，我们观察总体的速度差异，之后我们会观察估计中的误差。

>>> from sklearn.datasets import make_blobs 
>>> blobs, labels = make_blobs(int(1e6), 3)

>>> from sklearn.cluster import KMeans, MiniBatchKMeans

>>> kmeans = KMeans(n_clusters=3) >>> minibatch = MiniBatchKMeans(n_clusters=3)

要理解这些度量的目的是暴露问题。所以，需要多加小心，来确保跑分的高精度性。这个话题还有大量可用的信息。如果你真的希望了解，MiniBatch KMeans 为何在粒度上更具优势，最好还是要阅读它们。

既然准备已经完成，我们可以测量时间差异：

>>> %time kmeans.fit(blobs) #IPython Magic CPU times: user 8.17 s, sys: 881 ms, total: 9.05 s Wall time: 9.97 s

>>> %time minibatch.fit(blobs) CPU times: user 4.04 s, sys: 90.1 ms, total: 4.13 s Wall time: 4.69 s

CPU 时间上有很大差异。聚类性能上的差异在下面展示：

>>> kmeans.cluster_centers_[0] array([ 1.10522173, -5.59610761, -8.35565134])

>>> minibatch.cluster_centers_[0] array([ 1.12071187, -5.61215116, -8.32015587])

我们可能要问的下一个问题就是，两个形心距离多远。

>>> from sklearn.metrics import pairwise 
>>> pairwise.pairwise_distances(kmeans.cluster_centers_[0],
                                 minibatch.cluster_centers_[0])
array([[ 0.03305309]])

看起来十分接近了。对角线包含形心的差异：

>>> np.diag(pairwise.pairwise_distances(kmeans.cluster_centers_,
             minibatch.cluster_centers_)) 
array([ 0.04191979, 0.03133651,  0.04342707])

工作原理

这里的批次就是关键。批次被迭代来寻找批次均值。对于下一次迭代来说，前一个批次的均值根据当前迭代来更新。有多种选项，用于控制 KMeans 的通用行为，和决定 MiniBatch KMeans 的参数。

batch_size参数决定批次应为多大。只是玩玩的话，我们可以运行 MiniBatch，但是，此时我们将批次数量设置为和数据集大小相同。

>>> minibatch = MiniBatchKMeans(batch_size=len(blobs)) 
>>> %time minibatch.fit(blobs) 
CPU times: user 34.6 s, sys: 3.17 s, total: 37.8 s Wall time: 44.6 s

显然，这就违背了问题的核心，但是这的确展示了重要东西。选择差劲的初始条件可能影响我们的模型，特别是聚类模型的收敛。使用 MiniBatch KMeans，全局最优是否能达到，是不一定的。

3.5 使用 KMeans 聚类来量化图像

图像处理是个重要的话题，其中聚类有一些应用。值得指出的是，Python 中有几种非常不错的图像处理库。Scikit-image 是 Scikit-learn 的“姐妹”项目。如果你打算做任何复杂的事情，都值得看一看它。

准备

我们在这篇秘籍中会有一些乐趣。目标是使用聚类来把图像变模糊。

首先，我们要利用 SciPy 来读取图像。图像翻译为三维数组，x和y坐标描述了高度和宽度，第三个维度表示每个图像的 RGB 值。

# in your terminal 
$ wget http://blog.trenthauck.com/assets/headshot.jpg

操作步骤

现在，让我们在 Python 中读取图像：

>>> from scipy import ndimage 
>>> img = ndimage.imread("headshot.jpg") 
>>> plt.imshow(img)

下面就是图像：

嘿，这就是（年轻时期的）作者。

既然我们已经有了图像，让我们检查它的维度：

>>> img.shape 
(420, 420, 3)

为了实际量化图像，我们需要将其转换为二维数组，长为420x420，宽为 RGB 值。思考它的更好的方法，是拥有一堆三维空间中的数据点，并且对点进行聚类来降低图像中的不同颜色的数量 -- 这是一个简单的量化方式。

首先，让我们使数组变形，它是个 NumPy 数组，所以非常简单：

>>> x, y, z = img.shape 
>>> long_img = img.reshape(x*y, z) 
>>> long_img.shape (176400, 3)

现在我们开始聚类过程。首先，让我们导入聚类模块，并创建 KMeans 对象。我们传入n_clusters=5，使我们拥有 5 个簇，或者实际上是 5 个不同颜色。

这是个不错的秘籍，我们使用前面提到的轮廓距离：

>>> from sklearn import cluster 
>>> k_means = cluster.KMeans(n_clusters=5) 
>>> k_means.fit(long_img)

既然我们已经训练了 KMeans 对象，让我们看看我们的眼色：

>>> centers = k_means.cluster_centers_ 
>>> centers 
array([[ 142.58775848, 206.12712986,  226.04416873],
       [  86.29356543,  68.86312505,   54.04770507],
       [ 194.36182899,  172.19845258,  149.65603813],
       [  24.67768412,   20.45778933,   16.19698314],
       [ 149.27801776,  132.19850659,  115.32729167]])

工作原理

既然我们拥有了形心，我们需要的下一个东西就是标签。它会告诉我们，哪个点关联哪个簇。

>>> labels = k_means.labels_ 
>>> labels[:5] array([1, 1, 1, 1, 1], dtype=int32)

这个时候，我们需要最简的 NumPy 操作，之后是一个变形，我们就拥有的新的图像：

>>> plt.imshow(centers[labels].reshape(x, y, z))

下面就是产生的图像：

3.6 寻找特征空间中的最接近对象

有时，最简单的事情就是求出两个对象之间的距离。我们刚好需要寻找一些距离的度量，计算成对（Pairwise）距离，并将结果与我们的预期比较。

准备

Scikit-learn 中，有个叫做sklearn.metrics.pairwise的底层工具。它包含一些服务函数，计算矩阵X中向量之间的距离，或者X和Y中的向量距离。

这对于信息检索来说很实用。例如，提供一组客户信息，带有属性X，我们可能希望选取有个客户代表，并找到与这个客户最接近的客户。实际上，我们可能希望将客户按照相似性度量的概念，使用距离函数来排序。相似性的质量取决于特征空间选取，以及我们在空间上所做的任何变换。

操作步骤

我们会使用pairwise_distances函数来判断对象的接近程度。要记住，接近程度就像我们用于聚类/分类的距离函数。

首先，让我们从metric模块导入pairwise_distances函数，并创建用于操作的数据集：

>>> from sklearn.metrics import pairwise 
>>> from sklearn.datasets import make_blobs 
>>> points, labels = make_blobs()

用于检查距离的最简单方式是pairwise_distances：

>>> distances = pairwise.pairwise_distances(points)

distances是个 NxN的矩阵，对角线为 0。在最简单的情况中，让我们先看看每个点到第一个点的距离：

>>> np.diag(distances) [:5] 
array([ 0.,  0.,  0.,  0.,  0.])

现在我们可以查找最接近于第一个点的点：

>>> distances[0][:5] 
array([  0., 11.82643041,1.23751545, 1.17612135, 14.61927874])

将点按照接近程度排序，很容易使用np.argsort做到：

>>> ranks = np.argsort(distances[0]) 
>>> ranks[:5] 
array([ 0, 27, 98, 23, 67])

argsort的好处是，现在我们可以排序我们的points矩阵，来获得真实的点。

>>> points[ranks][:5] 
array([[ 8.96147382, -1.90405304],
       [ 8.75417014, -1.76289919],
       [ 8.78902665, -2.27859923],
       [ 8.59694131, -2.10057667],
       [ 8.70949958, -2.30040991]])

观察接近的点是什么样子，可能十分有用。结果在意料之中：

工作原理

给定一些距离函数，每个点都以成对函数来度量。通常为欧几里得距离，它是：

详细来说，它计算了两个向量每个分量的差，计算它们的平方，求和，之后计算它的平方根。这看起来很熟悉，因为在计算均方误差的时候，我们使用的东西很相似。如果我们计算了平方根，就一样了。实际上，经常使用的度量是均方根误差（RMSE），它就是距离函数的应用。

在 Python 中，这看起来是：

>>> def euclid_distances(x, y):
       return np.power(np.power(x - y, 2).sum(), .5) 
>>> euclid_distances(points[0], points[1])
11.826430406213145

Scikit-learn 中存在一些其他函数，但是 Scikit-learn 也会使用 SciPy 的距离函数。在本书编写之时，Scikit-learn 距离函数支持稀疏矩阵。距离函数的更多信息请查看 SciPy 文档。

cityblock
cosine
euclidean
l1
l2
manhattan

我们现在可以解决问题了。例如，如果我们站在原点处的格子上，并且线是街道，为了到达点(5,5)，我们需要走多远呢？

>>> pairwise.pairwise_distances([[0, 0], [5, 5]], metric='cityblock')[0] 
array([  0.,  10.])

使用成对距离，我们可以发现位向量之间的相似性。这是汉明距离的事情，它定义为：

使用下列命令：

>>> X = np.random.binomial(1, .5, size=(2, 4)).astype(np.bool) 
>>> X 
array([[False,  True, False, False],
       [False, False, False,  True]], dtype=bool)
>>> pairwise.pairwise_distances(X, metric='hamming') 
array([[ 0. ,  0.25],
       [ 0.25,  0. ]])

3.7 使用高斯混合模型的概率聚类

在 KMeans 中，我们假设簇的方差是相等的。这会导致空间的细分，这决定了簇如何被分配。但是，如果有一种场景，其中方差不是相等的，并且每个簇中的点拥有一个与之相关的概率，会怎么样？

准备

有一种更加概率化的方式，用于查看 KMeans 聚类。KMeans 聚类相当于将协方差矩阵S应用于高斯混合模型，这个矩阵可以分解为单位矩阵成误差。对于每个簇，协方差结构是相同的。这就产生了球形聚类。

但是，如果我们允许S变化，就可以估计 GMM，并将其用于预测。我们会以单变量的角度看到它的原理，之后扩展为多个维度。

操作步骤

首先，我们需要创建一些数据。例如，让我们模拟女性和男性的身高。我们会在整个秘籍中使用这个例子。这是个简单的例子，但是会展示出我们在 N 维空间中想要完成的东西，这比较易于可视化：

>>> import numpy as np 
>>> N = 1000

>>> in_m = 72 
>>> in_w = 66

>>> s_m = 2 
>>> s_w = s_m

>>> m = np.random.normal(in_m, s_m, N) 
>>> w = np.random.normal(in_w, s_w, N) 
>>> from matplotlib import pyplot as plt 
>>> f, ax = plt.subplots(figsize=(7, 5))

>>> ax.set_title("Histogram of Heights") 
>>> ax.hist(m, alpha=.5, label="Men"); 
>>> ax.hist(w, alpha=.5, label="Women"); 
>>> ax.legend()

下面是输出：

下面，我们的兴趣是，对分组二次抽样，训练分布，之后预测剩余分组。

>>> random_sample = np.random.choice([True, False], size=m.size) 
>>> m_test = m[random_sample] 
>>> m_train = m[~random_sample]

>>> w_test = w[random_sample] 
>>> w_train = w[~random_sample]

现在我们需要获得男性和女性高度的经验分布，基于训练集：

>>> from scipy import stats 
>>> m_pdf = stats.norm(m_train.mean(), m_train.std()) 
>>> w_pdf = stats.norm(w_train.mean(), w_train.std())

对于测试集，我们要计算，基于数据点从每个分布中生成的概率，并且最可能的分布会分配合适的标签。当然，我们会看到有多么准确。

>>> m_pdf.pdf(m[0]) 
0.043532673457165431
>>> w_pdf.pdf(m[0]) 
9.2341848872766183e-07

要注意概率中的差异。

假设当男性的概率更高时，我们会猜测，但是如果女性的概率更高，我们会覆盖它。

>>> guesses_m = np.ones_like(m_test) 
>>> guesses_m[m_pdf.pdf(m_test) < w_pdf.pdf(m_test)] = 0

显然，问题就是我们有多么准确。由于正确情况下guesses_m为 1，否则为 0，我们计算向量的均值来获取准确度。

>>> guesses_m.mean() 
0.93775100401606426

不是太糟。现在，来看看我们在女性的分组中做的有多好，使用下面的命令：

>>> guesses_w = np.ones_like(w_test) 
>>> guesses_w[m_pdf.pdf(w_test) > w_pdf.pdf(w_test)] = 0 
>>> guesses_w.mean() 0.93172690763052213

让我们允许两组间的方差不同。首先，创建一些新的数组：

>>> s_m = 1 
>>> s_w = 4

>>> m = np.random.normal(in_m, s_m, N) 
>>> w = np.random.normal(in_w, s_w, N)

之后，创建训练集：

>>> m_test = m[random_sample] 
>>> m_train = m[~random_sample]

>>> w_test = w[random_sample] 
>>> w_train = w[~random_sample] 
>>> f, ax = plt.subplots(figsize=(7, 5)) 
>>> ax.set_title("Histogram of Heights") 
>>> ax.hist(m_train, alpha=.5, label="Men"); 
>>> ax.hist(w_train, alpha=.5, label="Women"); 
>>> ax.legend()

让我们看看男性和女性之间的方差差异：

现在我们可以创建相同的 PDF：

>>> m_pdf = stats.norm(m_train.mean(), m_train.std()) 
>>> w_pdf = stats.norm(w_train.mean(), w_train.std())

下面是输出：

你可以在多维空间中想象他：

>>> class_A = np.random.normal(0, 1, size=(100, 2)) 
>>> class_B = np.random.normal(4, 1.5, size=(100, 2)) 
>>> f, ax = plt.subplots(figsize=(7, 5))

>>> ax.scatter(class_A[:,0], class_A[:,1], label='A', c='r')
>>> ax.scatter(class_B[:,0], class_B[:,1], label='B')

下面是输出：

工作原理

好的，所以既然我们看过了，我们基于分布对点分类的方式，让我们看看如何在 Scikit 中首先：

>>> from sklearn.mixture import GMM 
>>> gmm = GMM(n_components=2) 
>>> X = np.row_stack((class_A, class_B)) 
>>> y = np.hstack((np.ones(100), np.zeros(100)))

由于我们是小巧的数据科学家，我们创建训练集：

>>> train = np.random.choice([True, False], 200) 
>>> gmm.fit(X[train]) GMM(covariance_type='diag', init_params='wmc', min_covar=0.001,
     n_components=2, n_init=1, n_iter=100, params='wmc',
     random_state=None,   thresh=0.01)

训练和预测的完成方式，和 Scikit-learn 的其它对象相同。

>>> gmm.fit(X[train])
>>> gmm.predict(X[train])[:5] 
array([0, 0, 0, 0, 0])

既然模型已经训练了，有一些值得一看的其它方法。

例如，使用score_examples，我们实际上可以为每个标签获得每个样例的可能性。

3.8 将 KMeans 用于离群点检测

这一章中，我们会查看 Kmeans 离群点检测的机制和正义。它对于隔离一些类型的错误很实用，但是使用时应多加小心。

准备

这个秘籍中，我们会使用 KMeans，对簇中的点执行离群点检测。要注意，提及离群点和离群点检测时有很多“阵营”。以便面，我们可能通过移除离群点，来移除由数据生成过程生成的点。另一方面，离群点可能来源于测量误差或一些其它外部因素。

这就是争议的重点。这篇秘籍的剩余部分有关于寻找离群点。我们的假设是，我们移除离群点的选择是合理的。

离群点检测的操作是，查找簇的形心，之后通过点到形心的距离来识别潜在的离群点。

操作步骤

首先，我们会生成 100 个点的单个数据块，之后我们会识别 5 个离形心最远的点。它们就是潜在的离群点。

>>> from sklearn.datasets import make_blobs 
>>> X, labels = make_blobs(100, centers=1) 
>>> import numpy as np

非常重要的是，Kmeans 聚类只有一个形心。这个想法类似于用于离群点检测的单类 SVM。

>>> from sklearn.cluster import KMeans 
>>> kmeans = KMeans(n_clusters=1) 
>>> kmeans.fit(X)

现在，让我们观察绘图。对于那些远离中心的点，尝试猜测哪个点会识别为五个离群点之一：

>>> f, ax = plt.subplots(figsize=(7, 5)) 
>>> ax.set_title("Blob") 
>>> ax.scatter(X[:, 0], X[:, 1], label='Points') 
>>> ax.scatter(kmeans.cluster_centers_[:, 0],
                kmeans.cluster_centers_[:, 1],
                label='Centroid',
                color='r') 
>>> ax.legend()

下面就是输出：

现在，让我们识别五个最接近的点：

>>> distances = kmeans.transform(X) 
# argsort returns an array of indexes which will sort the array in ascending order 
# so we reverse it via [::-1] and take the top five with [:5] 
>>> sorted_idx = np.argsort(distances.ravel())[::-1][:5]

现在，让我们看看哪个点离得最远：

>>> f, ax = plt.subplots(figsize=(7, 5)) 
>>> ax.set_title("Single Cluster") 
>>> ax.scatter(X[:, 0], X[:, 1], label='Points') 
>>> ax.scatter(kmeans.cluster_centers_[:, 0],
                kmeans.cluster_centers_[:, 1],
                label='Centroid', color='r') 
>>> ax.scatter(X[sorted_idx][:, 0], X[sorted_idx][:, 1],
                label='Extreme Value', edgecolors='g',
                facecolors='none', s=100) 
>>> ax.legend(loc='best')

下面是输出：

如果我们喜欢的话，移除这些点很容易。

>>> new_X = np.delete(X, sorted_idx, axis=0)

同样，移除这些点之后，形心明显变化了。

>>> new_kmeans = KMeans(n_clusters=1) 
>>> new_kmeans.fit(new_X)

让我们将旧的和新的形心可视化：

>>> f, ax = plt.subplots(figsize=(7, 5)) 
>>> ax.set_title("Extreme Values Removed") 
>>> ax.scatter(new_X[:, 0], new_X[:, 1], label='Pruned Points') 
>>> ax.scatter(kmeans.cluster_centers_[:, 0],
               kmeans.cluster_centers_[:, 1], label='Old Centroid',
               color='r', s=80, alpha=.5) 
>>> ax.scatter(new_kmeans.cluster_centers_[:, 0],
               new_kmeans.cluster_centers_[:, 1], label='New Centroid',
               color='m', s=80, alpha=.5) 
>>> ax.legend(loc='best')

下面是输出：

显然，形心没有移动多少，仅仅移除五个极端点时，我们的预期就是这样。这个过程可以重复，知道我们对数据表示满意。

工作原理

我们已经看到，高斯分布和 KMeans 聚类之间有本质联系。让我们基于形心和样本的协方差矩阵创建一个经验高斯分布，并且查看每个点的概率 -- 理论上是我们溢出的五个点。这刚好展示了，我们实际上溢出了拥有最低可能性的值。距离和可能性之间的概念十分重要，并且在你的机器学习训练中会经常出现。

使用下列命令来创建经验高斯分布：

>>> from scipy import stats 
>>> emp_dist = stats.multivariate_normal(
               kmeans.cluster_centers_.ravel()) 
>>> lowest_prob_idx = np.argsort(emp_dist.pdf(X))[:5] 
>>> np.all(X[sorted_idx] == X[lowest_prob_idx]) True

3.9 将 KNN 用于回归

回归在这本书的其它地方有所设计，但是我们可能打算在特征空间的“口袋”中运行回归。我们可以认为，我们的数据集要经过多道数据处理工序。如果是这样，只训练相似数据点是个不错的想法。

准备

我们的老朋友，回归，可以用于聚类的上下文中。回归显然是个监督学习技巧，所以我们使用 KNN 而不是 KMeans。

对于 KNN 回归来说，我们使用特征空间中的 K 个最近点，来构建回归，而不像常规回归那样使用整个特征空间。

操作步骤

对于这个秘籍，我们使用iris数据集。如果我们打算预测一些东西，例如每朵花的花瓣宽度，根据iris物种来聚类可能会给我们更好的结果。KNN 回归不会根据物种来聚类，但我们的假设是，相同物种的 X 会接近，或者这个案例中，是花瓣长度。

对于这个秘籍，我们使用iris数据集：

>>> from sklearn import datasets 
>>> iris = datasets.load_iris() 
>>> iris.feature_names ['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)',  'petal width (cm)']

我们尝试基于萼片长度和宽度来预测花瓣长度。我们同时训练一个线性回归，来对比观察 KNN 回归有多好。

>>> from sklearn.linear_model import LinearRegression 
>>> lr = LinearRegression() 
>>> lr.fit(X, y) 
>>> print "The MSE is: {:.2}".format(np.power(y - lr.predict(X),
           2).mean()) 
The MSE is: 0.15

现在，对于 KNN 回归，使用下列代码：

>>> from sklearn.neighbors import KNeighborsRegressor 
>>> knnr = KNeighborsRegressor(n_neighbors=10) 
>>> knnr.fit(X, y) 
>>> print "The MSE is: {:.2}".format(np.power(y - knnr.predict(X),
           2).mean()) 
The MSE is: 0.069

让我们看看，当我们让它使用最接近的 10 个点用于回归时，KNN 回归会做什么？

>>> f, ax = plt.subplots(nrows=2, figsize=(7, 10))
>>> ax[0].set_title("Predictions")
>>> ax[0].scatter(X[:, 0], X[:, 1], s=lr.predict(X)*80, label='LR
    Predictions', color='c', edgecolors='black') 
>>> ax[1].scatter(X[:, 0], X[:, 1], s=knnr.predict(X)*80, label='k-NN
    Predictions', color='m', edgecolors='black')
>>> ax[0].legend() 
>>> ax[1].legend()

输出如下：

很显然，预测大部分都是接近的。但是让我们与实际情况相比，看看 Setosa 物种的预测：

>>> setosa_idx = np.where(iris.target_names=='setosa') 
>>> setosa_mask = iris.target == setosa_idx[0] 
>>> y[setosa_mask][:5] array([ 0.2,  0.2,  0.2,  0.2,  0.2]) 
>>> knnr.predict(X)[setosa_mask][:5] 
array([ 0.28,  0.17,  0.21,  0.2 ,  0.31])
>>> lr.predict(X)[setosa_mask][:5] 
array([ 0.44636645, 0.53893889, 0.29846368, 0.27338255, 0.32612885])

再次观察绘图，Setosa 物种（左上方的簇）被线性回归估计过高，但是 KNN 非常接近真实值。

工作原理

KNN 回归非常简单，它计算被测试点的 K 个最接近点的均值。

让我们手动预测单个点：

>>> example_point = X[0

现在，我们需要获取离我们的our_example_point最近的 10 个点：

>>> from sklearn.metrics import pairwise 
>>> distances_to_example = pairwise.pairwise_distances(X)[0] 
>>> ten_closest_points = X[np.argsort(distances_to_example)][:10] 
>>> ten_closest_y = y[np.argsort(distances_to_example)][:10]
>>> ten_closest_y.mean() 
0.28000

我们可以看到它非常接近预期。

Sklearn-秘籍-二-