编者按：机器学习开放课程第四课，Mail.Ru数据科学家Yury Kashnitsky深入讲解线性分类和线性回归的理论与实践。
xkcd译文：“如你所见，到下个月底，你会有48位丈夫。值得吃一大块蛋糕庆祝下。”（横轴为日期：昨天、今天；纵轴为丈夫数量。） 欢迎参加我们的第4课。这次我们将呈现最重要的主题——线性模型。如果你准备好了数据，打算开始训练模型，那么你最有可能首先尝试线性回归或逻辑回归（具体取决于你的任务是回归还是分类）。 本文包括线性模型的理论和实践（在实际任务中的使用）。此外，你将参加一项Kaggle竞赛，解决一个基于浏览历史识别用户的问题。 回归 最小二乘法 最大似然估计 偏置-方差分解 线性回归正则化
线性分类 线性分类器 作为线性分类器的逻辑回归 最大似然估计和逻辑回归 逻辑损失的L2正则化
逻辑回归正则化示例 逻辑回归的优缺点 验证和学习曲线 课内Kaggle竞赛“我知道你是谁” 作业四 相关资源
我们的线性模型学习从线性回归开始。首先，需要指定一个模型将因变量y和解释因素（特征）联系起来；对线性模型而言，依赖函数的形式如下： 
如果我们为每项观测加上一个虚维度x0 = 1 （偏置），那么上面的线性形式可以改写为一个略微紧凑的形式： 
如果我们有一个特征观测矩阵，其中矩阵的行是数据集中的观测，那么需要在左边加上一列。 我们可以定义模型为： 
其中： 上述表达式也可以写成这样（写出每项观察）： 
模型具有如下限制（否则它就不是线性回归了）： 权重wi的估计满足如下条件时，我们称其为线性： 
其中，∀k，ωki仅依赖于X中的样本，而且几乎一定以非线性的方式。由于寻求最佳权重的解是一个线性估计，这一模型被称为线性回归。让我们再引入一项定义。当期望值等于真实而未知的估计参数的值时，权重估计称为无偏（unbiased）： 
计算这些权重的方法之一是普通最小二乘法（OLS）。OLS最小化因变量的实际值和模型给出的预测值之间的均方误差： 
为了解决这一优化问题，我们需要计算模型参数的导数。我们将导数设为零，然后求解所得关于w的等式。 矩阵求导小抄： 
让我们开始计算： 
基于上述定义和条件，我们可以说，根据高斯-马尔可夫定理，模型参数的OLS估计是所有线性无偏估计中最优的，即，它们给出最低的方差。 最大似然估计 有人可能会问，为何我们选择最小化均方误差而不是别的什么？毕竟，我们可以最小化残差的平均绝对值。如果我们改变了最小化的值，唯一会发生的事就是我们将超出高斯-马尔可夫定理的条件，因此我们的估计将不再是最佳的线性无偏估计。 在我们继续之前，让我们稍微离题一下，通过一个简单的例子讲解下最大似然估计。 许多人大概都记得乙醇的化学式，所以我决定做一个试验判定人们是否记得简单的甲醇化学式：CH3OH. 我们调查了400人，发现只有117个人记得甲醇的化学式。那么，下一个受访者知道甲醇化学式的概率为117/400 ≈ 29%会是一个合理的假定。让我们展示下这个直观的估计不仅很好，同时也是最大似然估计。估计来自哪里？回忆下伯努利分布的定义：如果一个随机变量只有两个值（1和0，相应的概率为θ和1 - θ），那么该随机变量满足伯努利分布，遵循以下概率分布函数： 
这一分布正是我们所需要的，分布参数θ是一个人知道甲醇化学式的概率估计。在我们的400个独立试验中，让我们将试验的结果记为x = (x1, x2, ..., x400)。让我们写下数据（观测）的似然，即正好观测到117个随机变量θ = 1的实例和283个随机变量θ = 0的实例：
接着，我们将最大化这一θ的表达式。最常见的情况是，我们并不最大化似然p(x | θ)，转而最大化其对数（这一单调变换不影响解答，但大大简化了计算）：
为了找到最大化上式的θ，我们求θ的导数，设为零，求解所得等式：
结果发现，我们的直观估计正好是最大似然估计。现在让我们将这一推理过程应用到线性回归问题上，尝试找出均方误差背后有什么。为此，我们需要从概率论的角度来看线性回归。我们的模型和之前是一样的：
不过，现在让我们假定随机误差符合均值为零的正态分布：
据此改写模型：
由于样本是独立抽取的（不相关误差是高斯-马尔可夫定理的条件之一），数据的似然看起来会是密度函数p(yi)的积。让我们考虑对数似然，对数似然允许我们用和替换积：
我们想要找到最大似然假设，即，我们需要最大化表达式p(y ∣ X, w) 以得到wML。这和最大化其对数是一回事。注意，当我们针对某个参数最大化函数时，我们可以丢弃所有不依赖这一参数的成员：
所以，我们看到了，最大化数据的似然和最小化均方误差是一回事（给定以上的假定）。实际上，这是因为误差是正态分布的。 偏置-方差分解 让我们稍微谈下线性回归预测的误差性质（实际上，这一讨论在所有机器学习算法上都成立）。我们已经提到： 因此，点x的误差可分解为：
为了简明，我们将省略函数的参数。让我们分别考虑每个成员。据以下公式：
我们很容易就能分解前两项：
注意：
最后我们来处理和的最后一项。回忆一下，误差和目标变量相互独立：
最后，让我们把这些合并到一起：
我们已经达到了我们的终极目标——最后的等式告诉我们，任何线性模型的预测误差由三部分组成： 尽管我们对σ2无能为力，我们可以影响前两项。理想情况下，我们希望同时取消这两项（见下图中左上），但是，在实践中，常常需要在偏置和不稳定（高方差）间寻找平衡。
一般而言，当模型的计算增加了（例如，自由参数的数量增加了），估计的方差（分散程度）也会增加，但偏置会下降。由于模型完全记下了训练集而没能概括训练集，小小的变动将导致未预期的结果（过拟合）。另一方面，如果模型太弱，它将不能够学习模式，导致学习偏离正解较远的不同答案。
高斯-马尔可夫定理断言，在线性模型参数估计问题中，OLS估计是最佳的线性无偏估计。这意味着，如果存在任何无偏线性模型g，我们可以确信
线性回归正则化 在一些情形下，我们可能会为了稳定性（降低模型的方差）特意增加模型的偏置。高斯-马尔可夫定理的条件之一就是矩阵X是满秩的。否则，OLS解w = (XTX)-1XTy不存在，因为逆矩阵(XTX)-1不存在。换句话说，矩阵XTX将是奇异矩阵或退化矩阵。这被称为病态问题。这类问题必须加以矫正，也就是说，矩阵XTX需要变为非退化矩阵或非奇异矩阵（这正是这一过程叫做正则化的原因）。我们常常能在这类数据中观察到所谓的多重共线性：当两个或更多特征高度相关，也就是矩阵X的列之间“几乎”存在线性依赖。例如，在基于参数预测房价这一问题中，属性“含阳台面积”和“不含阳台面积”会有一个“几乎是”线性的关系。形式化地说，包含这类数据的矩阵XTX是可逆的，但由于多重共线性，一些本征值会接近零。在XTX的逆矩阵中，会出现一些极端巨大的本征值，因为逆矩阵的本征值为1/(λi)。这一本征值的波动会导致模型参数估计的不稳定，即，在训练数据中加入一组新的观测会导致完全不同的解。有一种正则化的方法称为吉洪诺夫正则化，大致上是在均方误差中加上一个新成员：
吉洪诺夫矩阵常常表达为单位矩阵乘上一个系数：
在这一情形下，最小化均方误差问题变为一个L2正则限定问题。如果我们对新的损失函数求导，设所得函数为零，据w重整等式，我们便得到了这一问题的解：
这类回归被称为岭回归。岭为对角矩阵，我们在XTX矩阵上加上这一对角矩阵，以确保我们能得到一个常规矩阵。
这样的解降低了分散程度，但增加了偏置，因为参数的正则向量同时最小化了，这导致解朝零移动。在下图中，OLS解为白色虚线的交点。蓝点表示岭回归的不同解。可以看到，通过增加正则化参数λ，我们使解朝零移动。
线性分类器 线性分类器背后的基本思路是，目标分类的值可以被特征空间中的一个超平面分开。如果这可以无误差地达成，那么训练集被称为线性可分。
我们之前已经了解了线性回归和普通最小二乘法（OLS）。现在考虑一个二元分类问题，将目标分类记为“+1”（正面样本）和“-1”（负面样本）。最简单的线性分类器可以通过回归定义：
其中 x是特征向量（包括标识）； w是线性模型中的权重向量（偏置为w0）； sign(∙)是符号函数，返回参数的符号； a(x)是分类x的分类器。
作为线性分类器的逻辑回归 逻辑回归是线性分类器的一个特殊情形，逻辑回归有一个额外的好处，可以预测样本xi为分类“+”的概率p+。
不仅能够预测一个回应（“+1”或“-1”），还能预测相应的概率，对于很多业务问题（比如，信用评分，这一问题传统上使用逻辑回归）而言，这是一个非常重要的需求。
银行选择一个阈值p*以预测贷款违约的概率（上图中阈值为0.15），超过阈值就不批准贷款。此外，还可以将预测概率乘以未偿还金额，以得到客户造成的期望损失，这可以构成有用的业务指标。 为了预测概率p+ ∈ [0, 1]，我们可以从使用OLS构造线性预测开始：
不过，为了将所得结果转换为[0, 1]区间内的概率，我们需要某个函数
逻辑回归使用如下函数：
%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np
def sigma(z):
return 1. / (1 + np.exp(-z))
xx = np.linspace(-10, 10, 1000)
plt.plot(xx, [sigma(x) for x in xx]);
plt.xlabel('z');
plt.ylabel('sigmoid(z)')
plt.title('Sigmoid function');
我们将事件X的概率记为P(X)，则比值比OR(X)由下式判定P(X)/(1-P(X))，这是某一事件是否发生的概率之比。显然，概率和比值比包含同样的信息，不过P(X)的范围是0到1，而OR(X)的范围是0到∞。 如果我们计算OR(X)的对数，那么显然我们有log OR(X) ∈ ℝ. 我们在OLS中将用到这个。 让我们看看逻辑回归是如何做出预测的：
目前而言，让我们假设我们已经通过某种方式得到了权重w，即，模型已经训练好了。之后我们将看下这是如何做的。 步骤一 计算
等式WTX = 0定义了将样本分为两类的超空间。 步骤二 计算对数比值比：
步骤三 现在我们已经有了将一个样本分配到“+”分类的概率OR+，我们可以据此计算p+：
我们看到，上式的右边我们得到了sigmoid函数。 所以，逻辑回归预测将一个样本分配为“+”分类的概率（假定我们已知模型的特征和权重），这一过程是通过对权重向量和特征向量的线性组合进行sigmoid变换完成的：
下面我们将看下模型是如何训练的。我们将再次依靠最大似然估计。 最大似然估计和逻辑回归 现在，让我们看下从MLE出发如何进行逻辑回归优化，也就是最小化逻辑损失函数。我们前面已经见过了将样本分配为“+”分类的逻辑回归模型：
“-”分类相应的表达式为：
这两个表达式可以组合成一个：
表达式M(xi) = yiwTxi称为目标xi的分类边缘。如果边缘非负，则模型正确选择了目标xi的分类；如果边缘为负，则目标xi被错误分类了。注意，边缘仅针对训练集中的目标（真实目标分类标签yi已知的目标）而言。 为了精确地理解我们为何得出这一结论，让我们转向线性分类器的几何解释。 首先，我会建议看下线性代数的一个经典入门问题：找出向径xA与平面wTx = 0的距离。 答案：
从答案中，我们可以看到，表达式wTxi的绝对值越大，点xi离平面wTx = 0的距离就越远。 因此，表达式M(xi) = yiwTxi是模型对目标xi分类的“信心”： 如果边缘的绝对值较大，且为正值，那么分类的标签是正确的，且目标离分界超平面很远，也就是模型对分类很自信。如下图点x3所示； 如果边缘的绝对值较大，且为负值，那么分类的标签是错误的，且目标离分界超平面很远（目标很可能是一个异常值；例如，它可能是训练集中一个错误标记的值）。如下图点x1所示； 如果边缘的绝对值较小，那么目标距离分界超平面很近，边缘的符号决定目标是否被正确分类了。如下图点x2和x4所示。
现在让我们计算数据集的似然，即基于数据集X观测到给定向量y的概率。我们将做一个强假设：目标来自一个独立分布（i.i.d.）。
其中，ℓ为数据集X的长度（行数）。 像我们经常干的那样，对这个表达式取对数，因为和要比积容易优化得多：
最大化似然等价于最小化以下表达式：
这就是逻辑损失函数。 用边缘改写逻辑损失函数，我们有：
我们将这一函数的图像和0-1损失函数的图像绘制在一张图上。0-1损失函数简单地惩罚模型误差（边缘为负）为1：
上图体现了这样一个想法，如果我们不能够直接最小化分类问题的误差数量（至少无法通过梯度方法最小化——0-1损失函数在零处的导数趋向无穷），我们可以最小化它的上界。对逻辑损失函数而言，以下是成立的：
我们希望能通过降低分类误差数的上界，降低分类误差数本身。 逻辑损失的L2正则化 逻辑回归的L2正则化和岭回归的情况基本一样。我们转而最小化下式：
在逻辑回归中，通常使用正则化系数的倒数C = 1/λ：
下面我们将通过一个例子直观地理解正则化。 让我们看下正则化是如何影响分类的质量的（数据集为吴恩达机器学习课程中的微芯片测试）。我们将使用基于多项式特征的逻辑回归，然后改变正则化参数C. 首先，我们将看看正则化是如何影响分类器的分界的，并直观地识别欠拟合和过拟合。接着，我们将通过交叉验证和网格搜索选择数值上接近最优值的正则化参数。 # 关闭警告
# 如果你喜欢开着警告，可以注释掉下面两行
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.model_selection import GridSearchCV
让我们使用pandas库的read_csv方法加载数据。这个数据集内有118个微芯片（目标），其中有两项质量控制测试的结果（两个数值变量）和微芯片是否投产的信息。变量已经居中了，也就是列中的值已经减去其均值。所以，“平均”微芯片的测试值为零。 # 加载数据
data = pd.read_csv('../../data/microchip_tests.txt',
header=None, names = ('test1','test2','released'))
# 了解数据集的基本信息
data.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 118 entries, 0 to 117
Data columns (total 3 columns):
test1 118 non-null float64
test2 118 non-null float64
released 118 non-null int64
dtypes: float64(2), int64(1)
memory usage: 2.8 KB
让我们看下开始五行和最后五行：
data.tail(5)
分离训练集和目标分类标签： X = data.iloc[:,:2].values
y = data.iloc[:,2].values
绘制数据，橙点对应有缺陷的芯片，蓝点对应正常芯片。 plt.scatter(X[y == 1, 0], X[y == 1, 1], c='blue', label='Released')
plt.scatter(X[y == 0, 0], X[y == 0, 1], c='orange', label='Faulty')
plt.xlabel('Test 1')
plt.ylabel('Test 2')
plt.title('2 tests of microchips. Logit with C=1')
plt.legend();
定义一个显示分类器的分界曲线的函数。 def plot_boundary(clf, X, y, grid_step=.01, poly_featurizer=None):
x_min, x_max = X[:, 0].min() - .1, X[:, 0].max() + .1
y_min, y_max = X[:, 1].min() - .1, X[:, 1].max() + .1
xx, yy = np.meshgrid(np.arange(x_min, x_max, grid_step),
np.arange(y_min, y_max, grid_step))
Z = clf.predict(poly_featurizer.transform(np.c_[xx.ravel(), yy.ravel()]))
Z = Z.reshape(xx.shape)
plt.contour(xx, yy, Z, cmap=plt.cm.Paired)
我们为两个变量x1和x2定义如下多形式特征：
例如，d = 3时的特征如下：
特征的数量呈指数型增长，为100个变量创建d较大（例如d = 10）的多项式特征成本很高。更重要的是，不需要如此。 我们将使用sklearn的逻辑回归实现。我们将创建一个对象，为矩阵X加上多项式特征（d不超过7）。 poly = PolynomialFeatures(degree=7)
X_poly = poly.fit_transform(X)
X_poly.shape
结果： 让我们训练逻辑回归，正则化系数C = 10-2。 C = 1e-2
logit = LogisticRegression(C=C, random_state=17)
logit.fit(X_poly, y)
plot_boundary(logit, X, y, grid_step=.01, poly_featurizer=poly)
plt.scatter(X[y == 1, 0], X[y == 1, 1], c='blue', label='Released')
plt.scatter(X[y == 0, 0], X[y == 0, 1], c='orange', label='Faulty')
plt.xlabel('Test 1')
plt.ylabel('Test 2')
plt.title('2 tests of microchips. Logit with C=%s' % C)
plt.legend();
print('Accuracy on training set:',
round(logit.score(X_poly, y), 3))
Accuracy on training set: 0.627
我们可以尝试将C增加到1，也就是说，我们削弱了正则化，现在的模型权重可以比之前有更大的值（绝对值更大）。这使得分类器在训练集上的精确度改善了（提高到0.831）。 C = 1
logit = LogisticRegression(C=C, random_state=17)
logit.fit(X_poly, y)
plot_boundary(logit, X, y, grid_step=.005, poly_featurizer=poly)
plt.scatter(X[y == 1, 0], X[y == 1, 1], c='blue', label='Released')
plt.scatter(X[y == 0, 0], X[y == 0, 1], c='orange', label='Faulty')
plt.xlabel('Test 1')
plt.ylabel('Test 2')
plt.title('2 tests of microchips. Logit with C=%s' % C)
plt.legend();
print('Accuracy on training set:',
round(logit.score(X_poly, y), 3))
Accuracy on training set: 0.831
我们为什么不接着进一步增加C呢？比如，将C增加到10000？这回，很明显正则化不够强，我们看到了过拟合。注意，C = 1时，“平滑”边界对应的训练集上的正确解答并没有低多少。但我们很容易可以想像得到，我们之前的模型将在新数据上工作得更好。 C = 1e4
logit = LogisticRegression(C=C, random_state=17)
logit.fit(X_poly, y)
plot_boundary(logit, X, y, grid_step=.005, poly_featurizer=poly)
plt.scatter(X[y == 1, 0], X[y == 1, 1], c='blue', label='Released')
plt.scatter(X[y == 0, 0], X[y == 0, 1], c='orange', label='Faulty')
plt.xlabel('Test 1')
plt.ylabel('Test 2')
plt.title('2 tests of microchips. Logit with C=%s' % C)
plt.legend();
print('Accuracy on training set:',
round(logit.score(X_poly, y), 3))
Accuracy on training set: 0.873
为了讨论以上这些结果，让我们改写一下逻辑回归优化的函数：
部分和： 参数C越大，模型可恢复的数据中的关系就越复杂（直观地说，C对应模型的“复杂度”——模型能力）。 如果正则化过强，即C值很小，最小化逻辑损失函数问题的解可能是一个许多权重过小或为零的解。这样的模型对误差的“惩罚”也不够（即，在函数J中，权重的平方和“权重过高”，误差L可能相对很大）。在这一情形下，模型将会欠拟合，如我们在第一个情形下所看到的那样。 相反，如果正则化过弱，即C值很大，由绝对值很大的分量组成的向量w可能变成优化问题的解。在这一情形下，L对优化的函数J贡献较大。大致上，这样的模型对在训练集的目标上犯错过于“恐惧”，因而会过拟合，如我们在第三个情形下所看到的那样。 逻辑回归不会“理解”（或“学习”）选择C的值，这和权重w的情况不一样。这就是说，C的值无法通过解决逻辑回归的优化问题而确定。我们之前碰到过类似的情况——决策树无法在训练过程中“学习”选择深度限制。因此，C是模型的超参数，通过交叉验证调节；决策树的max_depth同理。
正则化参数调整 让我们确定上述例子中正则化参数C的最优值。我们可以使用LogisticRegressionCV——网格搜索参数后进行交叉验证。这个类是专门为逻辑回归设计的。对任意模型而言，使用GridSearchCV或RandomizedSearchCV，或者特殊的超参数优化算法，比如hyperopt中实现的算法。 skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=17)
c_values = np.logspace(-2, 3, 500)
logit_searcher = LogisticRegressionCV(Cs=c_values, cv=skf, verbose=1, n_jobs=-1)
logit_searcher.fit(X_poly, y)
logit_searcher.C_
结果： 让我们看下超参数C是如何影响模型的质量的： plt.plot(c_values, np.mean(logit_searcher.scores_[1], axis=0))
plt.xlabel('C')
plt.ylabel('Mean CV-accuracy');
最后，选择C值“最佳”（C值较大加重算力负担，也容易导致过拟合）的区域： plt.plot(c_values, np.mean(logit_searcher.scores_[1], axis=0))
plt.xlabel('C')
plt.ylabel('Mean CV-accuracy');
plt.xlim((0,10));
回忆一下，这些曲线被称为验证曲线。之前我们手工创建了验证曲线，不过sklearn有构建这些曲线的特殊方法，我们以后将直接使用sklearn的相应方法。 分析IMDB影评 现在让我们做个小练习！我们想要解决IMDB影评的二元分类问题。我们有一个训练集，其中包含标记好的影评，12500条好评，12500条差评。直接开始机器学习并不容易，因为我们并没有矩阵X；我们需要准备它。我们将使用一个简单方法：词袋模型。影评的特征将由整个语料库中的每个词的出现情况表示。语料库是所有用户影评。下图展示了这一思路：
%matplotlib inline
import seaborn as sns
import numpy as np
from sklearn.datasets import load_files
from sklearn.feature_extraction.text import CountVectorizer, TfidfTransformer, TfidfVectorizer
from sklearn.linear_model import LogisticRegression
from sklearn.svm import LinearSVC
数据集可通过以下地址下载： http://ai./~amaas/data/sentiment/aclImdb_v1.tar.gz 数据集的简要描述： ai./~amaas/data/sentiment/ # 请修改文件路径
reviews_train = load_files('/Users/y.kashnitsky/Yandex.Disk.localized/ML/data/aclImdb/train')
text_train, y_train = reviews_train.data, reviews_train.target
看看训练集和测试集中各有多少条数据： print('Number of documents in training data: %d' % len(text_train))
print(np.bincount(y_train))
Number of documents in training data: 25000
[12500 12500]
# 请修改文件路径
reviews_test = load_files('/Users/y.kashnitsky/Yandex.Disk.localized/ML/data/aclImdb/test')
text_test, y_test = reviews_test.data, reviews_test.target
print('Number of documents in test data: %d' % len(text_test))
print(np.bincount(y_test))
Number of documents in test data: 25000
[12500 12500]
下面是影评的一些例子。 print(text_train[1])
b'Words can\'t describe how bad this movie is. I can\'t explain it by writing only. You have too see it for yourself to get at grip of how horrible a movie really can be. Not that I recommend you to do that. There are so many clich\xc3\xa9s, mistakes (and all other negative things you can imagine) here that will just make you cry. To start with the technical first, there are a LOT of mistakes regarding the airplane. I won\'t list them here, but just mention the coloring of the plane. They didn\'t even manage to show an airliner in the colors of a fictional airline, but instead used a 747 painted in the original Boeing livery. Very bad. The plot is stupid and has been done many times before, only much, much better. There are so many ridiculous moments here that i lost count of it really early. Also, I was on the bad guys\' side all the time in the movie, because the good guys were so stupid. 'Executive Decision' should without a doubt be you\'re choice over this one, even the 'Turbulence'-movies are better. In fact, every other movie in the world is better than this one.'
上面这条影评是差评还是好评？ y_train[1]
没错，和我们预料的一样，这是一条差评。 text_train[2]
b'Everyone plays their part pretty well in this 'little nice movie'. Belushi gets the chance to live part of his life differently, but ends up realizing that what he had was going to be just as good or maybe even better. The movie shows us that we ought to take advantage of the opportunities we have, not the ones we do not or cannot have. If U can get this movie on video for around $10, it\xc2\xb4d be an investment!'
这条呢？ y_train[2]
是好评。 单词简单计数 首先，我们使用CountVectorizer创建包含所有单词的字典。 cv = CountVectorizer()
cv.fit(text_train)
len(cv.vocabulary_)
如果你查看下“单词”的样本（我们还是称作token吧），你会发现我们省略了文本处理的许多重要步骤（自动化文本处理自身就可以成为一个独立的文章系列）。 print(cv.get_feature_names()[:50])
print(cv.get_feature_names()[50000:50050])
['00', '000', '0000000000001', '00001', '00015', '000s', '001', '003830', '006', '007', '0079', '0080', '0083', '0093638', '00am', '00pm', '00s', '01', '01pm', '02', '020410', '029', '03', '04', '041', '05', '050', '06', '06th', '07', '08', '087', '089', '08th', '09', '0f', '0ne', '0r', '0s', '10', '100', '1000', '1000000', '10000000000000', '1000lb', '1000s', '1001', '100b', '100k', '100m']
['pincher', 'pinchers', 'pinches', 'pinching', 'pinchot', 'pinciotti', 'pine', 'pineal', 'pineapple', 'pineapples', 'pines', 'pinet', 'pinetrees', 'pineyro', 'pinfall', 'pinfold', 'ping', 'pingo', 'pinhead', 'pinheads', 'pinho', 'pining', 'pinjar', 'pink', 'pinkerton', 'pinkett', 'pinkie', 'pinkins', 'pinkish', 'pinko', 'pinks', 'pinku', 'pinkus', 'pinky', 'pinnacle', 'pinnacles', 'pinned', 'pinning', 'pinnings', 'pinnochio', 'pinnocioesque', 'pino', 'pinocchio', 'pinochet', 'pinochets', 'pinoy', 'pinpoint', 'pinpoints', 'pins', 'pinsent']
接着，我们将使用单词的索引编码训练集文本的句子。我们将使用稀疏矩阵。 X_train = cv.transform(text_train)
X_train
<25000x74849 sparse matrix of type '<class 'numpy.int64'>'
with 3445861 stored elements in Compressed Sparse Row format>
让我们看下这一转换是如何进行的。 print(text_train[19726])
b'This movie is terrible but it has some good effects.'
X_train[19726].nonzero()[1]
array([ 9881, 21020, 28068, 29999, 34585, 34683, 44147, 61617, 66150, 66562], dtype=int32)
接下来，我们对测试集应用同样的操作。 X_test = cv.transform(text_test)
下一步是训练逻辑回归。 %%time
logit = LogisticRegression(n_jobs=-1, random_state=7)
logit.fit(X_train, y_train)
CPU times: user 40.9 s, sys: 524 ms, total: 41.4 s
Wall time: 10.5 s
让我们看下训练集和测试集上的精确度。 round(logit.score(X_train, y_train), 3), round(logit.score(X_test, y_test), 3)
(0.998, 0.86699999999999999)
模型的系数可以美观地显示。 def visualize_coefficients(classifier, feature_names, n_top_features=25):
# 获取绝对值较大的系数
coef = classifier.coef_.ravel()
positive_coefficients = np.argsort(coef)[-n_top_features:]
negative_coefficients = np.argsort(coef)[:n_top_features]
interesting_coefficients = np.hstack([negative_coefficients, positive_coefficients])
# 绘图
plt.figure(figsize=(15, 5))
colors = ['red' if c < 0 else 'blue' for c in coef[interesting_coefficients]]
plt.bar(np.arange(2 * n_top_features), coef[interesting_coefficients], color=colors)
feature_names = np.array(feature_names)
plt.xticks(np.arange(1, 1 + 2 * n_top_features), feature_names[interesting_coefficients], rotation=60, ha='right');
def plot_grid_scores(grid, param_name):
plt.plot(grid.param_grid[param_name], grid.cv_results_['mean_train_score'],
color='green', label='train')
plt.plot(grid.param_grid[param_name], grid.cv_results_['mean_test_score'],
color='red', label='test')
plt.legend();
visualize_coefficients(logit, cv.get_feature_names())
为了让我们的模型更好，我们可以优化逻辑回归的正则化系数。我们将使用sklearn.pipeline，因为CountVectorizer只应该应用于训练数据（为了避免“偷窥”测试集和在测试集上计算词频）。在这一情形下，pipeline确定正确的行动序列：应用CountVectorizer，然后训练逻辑回归。 %%time
from sklearn.pipeline import make_pipeline
text_pipe_logit = make_pipeline(CountVectorizer(),
LogisticRegression(n_jobs=-1, random_state=7))
text_pipe_logit.fit(text_train, y_train)
print(text_pipe_logit.score(text_test, y_test))
0.86672
CPU times: user 49.9 s, sys: 571 ms, total: 50.5 s
Wall time: 20.5 s
%%time
from sklearn.model_selection import GridSearchCV
param_grid_logit = {'logisticregression__C': np.logspace(-5, 0, 6)}
grid_logit = GridSearchCV(text_pipe_logit, param_grid_logit, cv=3, n_jobs=-1)
grid_logit.fit(text_train, y_train)
CPU times: user 26.7 s, sys: 1.33 s, total: 28.1 s
Wall time: 1min 16s
让我们查看下最佳C，以及相应的交叉验证评分： grid_logit.best_params_, grid_logit.best_score_
({'logisticregression__C': 0.10000000000000001}, 0.88527999999999996)
plot_grid_scores(grid_logit, 'logisticregression__C')
验证集上的结果： grid_logit.score(text_test, y_test)
0.87907999999999997
现在让我们用随机森林来分类。我们看到，逻辑回归事半功倍。 from sklearn.ensemble import RandomForestClassifier
forest = RandomForestClassifier(n_estimators=200, n_jobs=-1, random_state=17)
%%time
forest.fit(X_train, y_train)
CPU times: user 3min 39s, sys: 1.2 s, total: 3min 40s
Wall time: 30.7 s
RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
max_depth=None, max_features='auto', max_leaf_nodes=None,
min_impurity_split=1e-07, min_samples_leaf=1,
min_samples_split=2, min_weight_fraction_leaf=0.0,
n_estimators=200, n_jobs=-1, oob_score=False, random_state=17,
verbose=0, warm_start=False)
round(forest.score(X_test, y_test), 3)
XOR问题 现在让我们看一个线性模型表现不佳的例子。 线性分类定义的是一个非常简单的分界平面——一个超平面。最著名的超平面（或直线）无法无误差地切分的玩具例子是XOR问题。 XOR即异或，其真值表如下：
XOR是一个简单的二元分类问题，其中两个分类呈对角交叉分布。 # 创建数据集
rng = np.random.RandomState(0)
X = rng.randn(200, 2)
y = np.logical_xor(X[:, 0] > 0, X[:, 1] > 0)
plt.scatter(X[:, 0], X[:, 1], s=30, c=y, cmap=plt.cm.Paired);
显然，我们无法划出一条直线无误差地将两个分类分开。因此，逻辑回归在这一任务上的表现很差。 def plot_boundary(clf, X, y, plot_title):
xx, yy = np.meshgrid(np.linspace(-3, 3, 50),
np.linspace(-3, 3, 50))
clf.fit(X, y)
# 为网格上的每个数据点绘制判定函数
Z = clf.predict_proba(np.vstack((xx.ravel(), yy.ravel())).T)[:, 1]
Z = Z.reshape(xx.shape)
image = plt.imshow(Z, interpolation='nearest',
extent=(xx.min(), xx.max(), yy.min(), yy.max()),
aspect='auto', origin='lower', cmap=plt.cm.PuOr_r)
contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,
linetypes='--')
plt.scatter(X[:, 0], X[:, 1], s=30, c=y, cmap=plt.cm.Paired)
plt.xticks(())
plt.yticks(())
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.axis([-3, 3, -3, 3])
plt.colorbar(image)
plt.title(plot_title, fontsize=12);
plot_boundary(LogisticRegression(), X, y,
'Logistic Regression, XOR problem')
然而，如果我们给出多项式特征作为输入（这里d = 2），那么问题解决了。 from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
logit_pipe = Pipeline([('poly', PolynomialFeatures(degree=2)),
('logit', LogisticRegression())])
plot_boundary(logit_pipe, X, y,
'Logistic Regression + quadratic features. XOR problem')
这里，逻辑回归仍然生成了一个超平面，不过是6维特征空间（1、x1、x2、x12、x1x2、x22）中的超平面。当我们将这个超平面投影到原特征空间（x1、x2）时，分界是非线性的。 在实践中，多项式特征确实有帮助，不过显式创建它们在算力上是低效的。使用核技巧的SVM要快很多。在这一方法中，只计算高维空间中目标之间的距离（由核函数定义），而不用生成大量特征组合。 现在，我们对模型验证、交叉验证、正则化已经有所了解，让我们考虑一个更大的问题： 如果模型的质量不尽人意，该怎么办？ 我们应该让模型更复杂还是更简单？ 我们应该加入更多特征吗？ 我们是否只是需要更多数据用于训练？
这些问题的答案并不明显。特别是，有时候一个更复杂的模型会导致表现退化。另一些时候，增加新的观测并不会带来可以观察得到的变化。事实上，做出正确决定，选择正确方法，从而改进模型的能力区分了优秀的专业人员和糟糕的专业人员。 让我们回头看看电信运营商离网率数据集。 %matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV, SGDClassifier
from sklearn.model_selection import validation_curve
data = pd.read_csv('../../data/telecom_churn.csv').drop('State', axis=1)
data['International plan'] = data['International plan'].map({'Yes': 1, 'No': 0})
data['Voice mail plan'] = data['Voice mail plan'].map({'Yes': 1, 'No': 0})
y = data['Churn'].astype('int').values
X = data.drop('Churn', axis=1).values
我们将使用随机梯度下降训练逻辑回归。我们将在之后的课程专门讨论梯度下降。 alphas = np.logspace(-2, 0, 20)
sgd_logit = SGDClassifier(loss='log', n_jobs=-1, random_state=17)
logit_pipe = Pipeline([('scaler', StandardScaler()), ('poly', PolynomialFeatures(degree=2)),
('sgd_logit', sgd_logit)])
val_train, val_test = validation_curve(logit_pipe, X, y,
'sgd_logit__alpha', alphas, cv=5,
scoring='roc_auc')
我们将绘制ROC-AUC曲线，查看下不同正则化参数导致的模型在训练集和参数集上的不同表现。
趋势相当明显，这样的趋势在其他问题中也同样常见： 需要多少数据？ 数据是越多越好。但我们如何理解新数据是否在任何情况下都有帮助呢？例如，如何评估花费N来加倍数据集是否合理呢？ 由于新数据可能无法取得，合理的做法是改变训练集的大小，然后看解答的质量如何依赖于训练数据的数量。这样我们就得到了学习曲线。 这个想法很简单：我们将误差显示为训练中使用的样本数的函数。模型的参数事先固定。 from sklearn.model_selection import learning_curve
def plot_learning_curve(degree=2, alpha=0.01):
train_sizes = np.linspace(0.05, 1, 20)
logit_pipe = Pipeline([('scaler', StandardScaler()), ('poly', PolynomialFeatures(degree=degree)),
('sgd_logit', SGDClassifier(n_jobs=-1, random_state=17, alpha=alpha))])
N_train, val_train, val_test = learning_curve(logit_pipe,
X, y, train_sizes=train_sizes, cv=5,
scoring='roc_auc')
plot_with_err(N_train, val_train, label='training scores')
plot_with_err(N_train, val_test, label='validation scores')
plt.xlabel('Training Set Size'); plt.ylabel('AUC')
plt.legend()
让我们看看线性模型的结果。我们将把正则化系数设定为较大的数。 plot_learning_curve(degree=2, alpha=10)
一个典型的状况：对于少量数据而言，训练集和交叉验证集之间的误差差别相当大，这暗示了过拟合。同样的模型，不过使用大量数据，误差“收敛”，暗示了欠拟合。 如果我们加入更多数据，训练集的误差不会增加。另一方面，测试集上的误差也不会下降。 所以，我们看到误差“收敛”，加入新数据无济于事。事实上这个情况对业务来说很重要。我们可能把数据集大小增大10倍，但是，如果我们不改变模型的复杂度，数据的增加没有帮助。因此“设定一次，使用十次”的策略可能无法奏效。 如果我们将正则化系数降低到0.05，会怎么样？ 我们看到了好兆头——曲线逐渐收敛，如果我们移动到更右，也就是加入更多数据，我们可以进一步改善模型在验证集上的表现。 plot_learning_curve(degree=2, alpha=0.05)
如果我们把alpha设为10-4，让模型更复杂，会出现什么情况？ 我们看到了过拟合——在训练集和验证集上，AUC都下降了。 plot_learning_curve(degree=2, alpha=1e-4)
构建这些曲线可以帮助我们理解朝哪个方向走，如何恰当地为新数据调整模型复杂度。 学习曲线和验证曲线的一些结论： 训练集上的误差本身不能说明模型的质量。 交叉验证误差显示了模型对数据拟合得有多好（数据的现有趋势）同时保留了多少对新数据的概括能力。 验证曲线是一个根据模型复杂度显示训练集和验证集上的结果的图形： 学习曲线是一个根据观测数量显示训练集和验证集上的结果的图形：
Kaggle竞赛页面： https://www./c/catch-me-if-you-can-intruder-detection-through-webpage-session-tracking2 我们将解决一个入侵者检测问题（通过分析上网用户的行为）。这是一个结合了数据分析和行为心理学的复杂而有趣的问题。例如，Yandex根据用户的行为模式解决邮箱入侵检测问题。概括起来，入侵者的行为模式可能和邮箱所有者不一样： 所以我们可以检测到入侵者，将其扔出邮箱，要求通过短信验证码认证身份。 Google Analytics研发了类似的技术，并发表了相应的论文。你可以通过搜索“Traversal Pattern Mining”（遍历模式挖掘）和“Sequential Pattern Mining”（序列模式挖掘）了解更多关于这一主题的信息。 在这一竞赛中，我们将解决一个类似的问题：我们的算法将分析特定用户频繁访问的网站序列，预测这个人是不是一个名为Alice的用户还是一个入侵者（其他人）。我们将使用ROC-AUC作为指标。 数据下载和转换 如果你之前没有注册过Kaggle账号的话，注册一下。然后到竞赛页面，下载数据。 首先，加载训练集和测试集。探索下数据： %matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
import pickle
import numpy as np
import pandas as pd
from scipy.sparse import csr_matrix
from scipy.sparse import hstack
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import roc_auc_score
from sklearn.linear_model import LogisticRegression
# 加载训练集和测试集
train_df = pd.read_csv('../../data/websites_train_sessions.csv',
index_col='session_id')
test_df = pd.read_csv('../../data/websites_test_sessions.csv',
index_col='session_id')
# 转换time1、……、time10列至datetime类
times = ['time%s' % i for i in range(1, 11)]
train_df[times] = train_df[times].apply(pd.to_datetime)
test_df[times] = test_df[times].apply(pd.to_datetime)
# 据时间排序数据
train_df = train_df.sort_values(by='time1')
# 查看下训练集的开始几行
train_df.head()
训练集包含以下特征： 用户会话的选取方式保证这些会话不超过半小时或包含超过十个网站，也就是说，一旦一个用户访问过了十个网站，或者会话时长超过了半小时，我们就认为这一会话结束了。 表中有一些空值，意味着某些会话包含不到十个网站。将这些空值替换为0，并将列类型改为整数。同时，加载网站字典，看看是什么样的： # 将site1、……、site10列类型转为整数，同时用零填充NA值
sites = ['site%s' % i for i in range(1, 11)]
train_df[sites] = train_df[sites].fillna(0).astype('int')
test_df[sites] = test_df[sites].fillna(0).astype('int')
# 加载网站字典
with open(r'../../data/site_dic.pkl', 'rb') as input_file:
site_dict = pickle.load(input_file)
# 创建字典的dataframe
sites_dict = pd.DataFrame(list(site_dict.keys()),
index=list(site_dict.values()), columns=['site'])
print(u'Websites total:', sites_dict.shape[0])
sites_dict.head()
简单起见，我们将只使用会话中访问的站点（不考虑访问时长）。这一数据选择背后的依据是：Alice有她偏爱的站点，我们在会话中看到更多这些站点，这一会话是Alice的会话的可能性就越高，反之亦然。 让我们准备下数据。首先，我们从训练集中排除目标变量，这样，训练集和测试集就有了相同数量的列，我们可以将其合并为一个dataframe，对其一并进行转换。 # 目标变量
y_train = train_df['target']
# 合并为一个dataframe
full_df = pd.concat([train_df.drop('target', axis=1), test_df])
# 分割训练集和测试集的索引
idx_split = train_df.shape[0]
只保留dataframe的site1, site2, ..., site10。 # 访问过的网站索引构成的dataframe
full_sites = full_df[sites]
full_sites.head()
会话是网站索引的序列，这一数据表示不便于通过线性方法处理。我们需要将其转换为如下表示形式：每个网站对应一个特征（列），其值为会话中访问该网站的次数。 # 索引序列
sites_flatten = full_sites.values.flatten()
# 我们想要的矩阵
full_sites_sparse = csr_matrix(([1] * sites_flatten.shape[0],
sites_flatten,
range(0, sites_flatten.shape[0] + 10, 10)))[:, 1:]
如果你明白刚刚发生了什么，那你可以直接跳到下一节（也许你也可以处理逻辑回归？），否则，让我们一起来搞明白发生了什么。 稀疏矩阵 让我们估计一下，在上面的例子中，储存我们的数据需要多少内存。合并后的dataframe包含336k样本，每个样本包含48k整数特征。因此我们大致需要的内存为： 336K * 48K * 8 bytes = 16G * 8 Bytes = 128 GB
显然，凡夫俗子没有这么多内存（严格来说，Python可能允许你创建这样一个矩阵，但对其进行任何操作都不容易）。一个有趣的事实是，我们的矩阵的大多数元素值为零。如果我们计算非零元素，那么大概有一百八十万，所占比例略高于所有矩阵元素的0.01%. 这样一种大多数元素为零的矩阵，称为稀疏矩阵，零元素数量和总元素的比例则称为矩阵的稀疏度。 scipy.sparse库可用于处理稀疏矩阵，参考它的文档了解稀疏矩阵的类型，如何处理，何时使用稀疏矩阵最高效。注意，稀疏矩阵只包含非零元素。最后，让我们看看稀疏矩阵占用的内存大小（显然，稀疏矩阵显著节省了内存）： # 稀疏矩阵占用多少内存？
print('{0} elements * {1} bytes = {2} bytes'.
format(full_sites_sparse.count_nonzero(), 8,
full_sites_sparse.count_nonzero() * 8))
# 或者直接：
print('sparse_matrix_size = {0} bytes'.
format(full_sites_sparse.data.nbytes))
1866898 elements * 8 bytes = 14935184 bytes
sparse_matrix_size = 14935184 bytes
让我们通过一个迷你的例子探索下这一矩阵是如何形成的。假设我们的用户会话表是这样的：
总共有3个会话，每次会话不超过3个站点。用户访问的不同站点总数为4（表中的1到4表示这4个站点）。让我们假定这4个站点是： vk.com habrahabr.ru yandex.ru ods.ai
如果用户在会话时访问了不到三个站点，后面的值将为零。我们想要将原dataframe转换为每个会话显示某一特定站点的访问数，如下表所示：
为此，我们使用csr_matrix ((data, indices, indptr))创建一个频率表。这里，我们显式地设定所有参数，让你能更清楚地理解这一过程： # 创建由1组成的列表，长度等于原dataframe中元素的数量（9）
data = [1] * 9
# 网站id
indices = [1, 0, 0, 1, 3, 1, 2, 3, 4]
# 切分行（会话）的索引
indptr = [0, 3, 6, 9]
# 将上述三个变量聚合为一个元组，然后构建矩阵
# 显示矩阵时，将其转为通常的“密集”矩阵
csr_matrix((data, indices, indptr)).todense()
matrix([[2, 1, 0, 0, 0],
[0, 2, 0, 1, 0],
[0, 0, 1, 1, 1]])
你也许注意到了，所得矩阵的列数目不是四（不同网站的总数），而是五。第零列是额外增加的列，显示每次会话少访问了几个站点（在我们的迷你例子中，会话的长度为3）。这一列是多余的，需要从dataframe中移除。 使用稀疏矩阵的另一个好处是，有为其特制的矩阵操作和机器学习算法实现，有时能通过利用数据结构的特点显著加速操作。逻辑回归也适用。现在，万事俱备，我们可以创建第一个模型了。 训练第一个模型 我们将使用sklearn的逻辑回归实现（默认参数）。数据集的前90%用于训练（数据集按时间排序），剩余10%用于验证。 def get_auc_lr_valid(X, y, C=1.0, seed=17, ratio = 0.9):
# 将数据分为训练集和验证集
idx = int(round(X.shape[0] * ratio))
# 训练分类器
lr = LogisticRegression(C=C, random_state=seed,
solver='lbfgs', n_jobs=-1).fit(X[:idx, :], y[:idx])
# 为验证集做出预测
y_pred = lr.predict_proba(X[idx:, :])[:, 1]
# 计算质量
score = roc_auc_score(y[idx:], y_pred)
return score
%%time
# 选择训练集
X_train = full_sites_sparse[:idx_split, :]
# 在验证集上计算量度
print(get_auc_lr_valid(X_train, y_train))
0.9198622553850315
CPU times: user 138 ms, sys: 77.1 ms, total: 216 ms
Wall time: 2.74 s
第一个模型在验证集上的表现接近0.92 ROC-AUC。让我们将其作为第一条基线（作为一个开始）。为了在测试集上进行预测，我们需要在整个训练集上再次训练模型（到此为止，我们的模型只使用了数据的一部分进行训练），这将增加模型的概括能力： # 将预测写入文件的函数
def write_to_submission_file(predicted_labels, out_file,
target='target', index_label='session_id'):
predicted_df = pd.DataFrame(predicted_labels,
index = np.arange(1,
predicted_labels.shape[0] + 1),
columns=[target])
predicted_df.to_csv(out_file, index_label=index_label)
# 在整个训练数据集上训练模型
# 为了可重复性，将random_state设为17
# 显式设置C=1（这是默认值）
lr = LogisticRegression(C=1.0, solver='lbfgs',
random_state=17).fit(X_train, y_train)
# 在测试数据集上进行预测
X_test = full_sites_sparse[idx_split:,:]
y_test = lr.predict_proba(X_test)[:, 1]
# 写入预测结果至提交文件
write_to_submission_file(y_test, 'baseline_1.csv')
如果你遵循这些步骤，并将答案上传到竞赛页面，你将在公开排行榜上取得ROC AUC = 0.91707的成绩。 你的任务是通过特征工程、特征缩放和正则化进一步改善模型。你将首先尝试添加一些明显的特征（比如浏览网站的时刻，会话中的网站数目，等等）。我们鼓励你在课程的学习过程中尝试新的想法和模型，并参与竞赛——这很有趣！ 截止日期：March 11, 23:59 CET
I. Goodfellow、Y. Bengio、A. Courville所著《Deep Learning》（深度学习）一书提供了紧凑而出色的线性模型综述。 基本上每本ML教材都涉及线性模型。我们推荐C. Bishop的《Pattern Recognition and Machine Learning》和K. Murphy的《Machine Learning: A Probabilistic Perspective》。 如果你打算从统计学的视角概览线性模型，可以看下T. Hastie、R. Tibshirani、J. Friedman的《The elements of statistical learning》。 P. Harrington的《Machine Learning in Action》将引导你完全使用Python实现经典ML算法。 scikit-learn库。scikit-learn的开发者致力于编写极为清晰的文档。 Scipy 2017 scikit-learn教程（Alex Gramfort、Andreas Mueller）。 MTH594课程 Advanced data mining: theory and applications包含很多非常好的材料。 GitHub仓库rushter/MLAlgorithms里有许多ML算法的实现，其中包括线性回归和逻辑回归。 欢迎留言分享其他资源。
原文地址：https:///open-machine-learning-course/open-machine-learning-course-topic-4-linear-classification-and-regression-44a41b9b5220
|
