MNIST 数据集 (Mixed National Institute of Standards and Technology database) 是美国国家标准与技术研究院收集整理的大型手写数字数据库,包含60,000个示例的训练集以及10,000个示例的测试集。

获取数据

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

    memory = joblib.Memory('./view')
    fetch_openml_cached = memory.cache(fetch_openml)
    mnist_dataset = fetch_openml_cached('mnist_784', as_frame=True)
    mnist_dataset = pd.DataFrame(data=np.c_[mnist_dataset['data'], mnist_dataset['target']],
                                 columns=mnist_dataset['feature_names'] + ['target'])
    return mnist_dataset

打包为函数 load_mnist

查看分析

pd.info()

pd.describe()项的含义和相应的单个方法

• count：元素数
• unique：具有唯一（唯一）值的元素数
• top： 最频繁值（mode）
• freq：频率（出现次数）
• mean： 算术平均值
• std： 标准差
• min：最小值
• max： 最大值
• 50%： 中位数（median）
• 25%， ： 1/4 分位数， 3/4 分位数75%

    print(mnist_dataset)
    mnist_dataset.info()

    print(mnist_dataset.columns)
    print(mnist_dataset.describe())

数据预处理

自己定义查看 numpy 数据

def view_data(X, line_num=10):
    print("\nview_data\nshape:{0}, dim: {1}, dtype: {2}".format(X.shape, X.ndim, X.dtype))
    print(X[:line_num])

包装函数

    X, y = mnist_matrix[:, 0:-1], mnist_matrix[:, -1]
    view_data(X, 1)
    view_data(y)

转换为合适的数据结构

    mnist_matrix = np.array(mnist_matrix, dtype=float)

查看数据、可视化

    digit = np.array(X[20], dtype=int)
    digit = digit.reshape(28, 28)
    plt.imshow(digit, cmap="binary")
    plt.show()

    plt.hist(y)
    plt.hist(y[:60000])
    plt.show()

归一化

    X = mnist_matrix[:, 0:-1] / 255
    y = np.array(mnist_matrix[:, -1], dtype=int)
    view_data(X, 1)
    view_data(y)

划分数据集

    train_threshold = 60000
    N = y.shape[0]
    X = np.hstack((np.ones((N, 1)), X))

    X_training = X[:train_threshold, :]
    y_training = y[:train_threshold]

    X_test = X[train_threshold:, :]
    y_test = y[train_threshold:]

训练

获取维度、进一步处理训练集、初始化模型、误差函数

批量损失函数是这样的：

\mathrm{Lost}\left(h_{\theta}(x), y\right)=\left\{\begin{array}{ll}
-\log \left(h_{\theta}(x)\right) & y=1 \\
-\log \left(1-h_{\theta}(x)\right) & y=0
\end{array}\right.

合并为这个式子：

\mathrm{Lost}\left(h_{\theta}(x), y\right)=-y \log \left(h_{\theta}(x)\right)-(1-y) \log \left(1-h_{\theta}(x)\right)

模型

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

def hypothesis(W, X):
    return sigmoid(X @ W)

def logistic_regression(X, y, aim_digit):
    N, dim = X.shape[0], X.shape[1]
    W = np.zeros(dim, dtype=np.float32)

    tmp = aim_digit * np.ones(N, dtype=np.int32)
    y_new = np.array((y == tmp), dtype=np.int32)
    # view_data(y)

    error_func(W, X, y_new)
    result = optimize.minimize(fun=error_func, x0=W, args=(X, y), jac=gradient,  method='Newton-CG')
    return result.x

测试、评价模型

试一下

    count = 5
    X_sample = X_test[:count, :]
    for i in range(count):
        digit = np.array(X_sample[i, 1:] * 255, dtype=int)
        digit = digit.reshape(28, 28)
        plt.imshow(digit, cmap="binary")
        plt.show()
    y_hat = predict(W, X_sample)
    print(y_hat == 1)

评价指标

准确率(accuracy)计算公式为：

A C C=\frac{T P+T N}{T P+T N+F P+F N}

精确率是针对我们预测结果而言的，它表示的是预测为正的样本中有多少是真正的正样本。你认为的正样本，有多少猜对了（猜的精确性如何）。那么预测为正就有两种可能了，一种就是把正类预测为正类(TP)，另一种就是把负类预测为正类(FP)，也就是

P=\frac{TP}{TP+FP}

而召回率表示的是样本中的正例有多少被预测正确了。或者说正样本有多少被找出来了（召回了多少）。那也有两种可能，一种是把原来的正类预测成正类(TP)，另一种就是把原来的正类预测为负类(FN)。

R=\frac{TP}{TP+FN}

在信息检索领域，精确率和召回率又被称为查准率和查全率，

查准率＝检索出的相关信息量 / 检索出的信息总量 \\
查全率＝检索出的相关信息量 / 系统中的相关信息总量

F-Measure 是 Precision 和 Recall 加权调和平均。F1 度量为

F = \frac{2\times P\times R}{P+R}

def get_accuracy(W, X_test, y_test, aim_digit):
    N = y_test.shape[0]
    y_test = y_test == aim_digit
    y_hat = hypothesis(W, X_test) > 0.5
    view_data(y_test)
    view_data(y_hat)

    tmp = np.array(y_test == y_hat, dtype=np.int32)
    view_data(tmp)
    return np.mean(tmp) * 100

建立多分类器

利用 sigmoid 函数输出各个不同类型数字的概率，取最高概率，就得到了一个预测数字。有兴趣可以实现一下。

scikit-learn

from sklearn.linear_model import LogisticRegression
logisticRegr = LogisticRegression(solver="newton-cg")
logisticRegr.fit(X_training, y_training)
y_hat = logisticRegr.predict(X_test)
score = logisticRegr.score(X_test, y_test)
print(score)

结构

• 决策树原理

  • 决策理论
  • 学习过程

• 特征选择与决策树生成

  • 熵模型：信息与信息熵
  • ID3 算法：信息增益
  • C4.5 算法：信息增益比
  • CART 算法：基尼系数
  • 决策树剪枝

决策树原理

来自人的决策经验

决策树类似我们人的决策过程。比如买北京学区房。就读学校是什么？预算？面积？

x = {x_1, x_2, \dots,x_n}

决策树算法利用树形结构，使用层层递进提问得到最终的分类。通过上面的问题，我们就构建了一套买房的策略。

• 每个非叶节点表示一个特征属性测试。
• 每个分支代表这个特征属性在某个值域（分类）上的输出。
• 每个叶子节点存放一个类别。
• 每个节点包含的样本集合通过属性测试被划分到子节点中，根节点包含样本全集。

学习过程

决策树包含根节点、内部节点和叶节点，其根节点包含样本全集，内部节点对应特征属性测试，叶节点则代表决策结果。

由于决策树是基于特征对实例进行分类的，因而其学习的本质是从训练数据集中归纳出一组用于分类的 if-else 规则。在学习的过程中，这组规则集合既要在训练数据上降低经验误差，也要具备良好的泛化能力（降低泛化误差）。

决策树模型的学习过程包括三个步骤：特征选择、决策树生成、决策树剪枝。

特征选择与决策树生成

熵模型：信息与信息熵

特征选择决定了使用哪些特征来划分特征空间。训练集中，每个样本的属性可能有很多个，在分类结果中起到的作用也有大有小。因此特征选择作用在于筛选出与分类结果相关性较高，也就是分类能力较强的特征。理想的特征选择是在每次划分之后，分支节点所包含的样本都尽可能属于同一个类别。

在特征选择中通常使用的准则是信息增益。信息增益描述的是在已知特征后对数据分类不确定性的减少程度，因而特征的信息增益越大，得到的分类结果的不确定度越低，特征也就具有越强的分类能力。根据信息增益准则选择特征的过程，就是自顶向下进行划分，在每次划分时计算每个特征的信息增益并选取最大值的过程。信息增益的计算涉及信源熵和条件熵的公式。

如果事件 X 发生的概率是 p(X)，这个事件的自信息量的定义为

\mathrm{H}(X) = -log_2(p(X))

什么是信息，如何度量信息？

选择题有选项 A/B/C/D，结合我们朴素的生活观念，以下那句话带给我们的信息量最大？

1. 「D一定是错误答案。」
2. 「每个选项的概率是均等的。」
3. 「正确答案是 A 。」

可以知道越大概率的事情发生了产生的信息量越小。

参考 YJango 的视频

度量质量需要一个物体的质量作为单位，自然地，我们也想到用一件事（随机变量）的不确定性作为另一件事的单位。

比如一枚硬币的可能性是 2，也就是两种不确定性，每种情况发生的概率是 1/2。

为了可比性，3 枚硬币的不确定性应该是一枚硬币不确定性的三倍，而 3 枚硬币的不确定性是 8。为了方便我们比较，我们选择用对数运算。

log_2(2^3) = 3

所以 3 枚硬币的不确定性，若以一枚硬币的不确定性（当选择的参照物有两个等概率事件）作为单位，它的信息熵是 3 bit 。

那么选择题的 ABCD 选项在我们不知道任何提示的情况下（ ABCD 是四种等可能事件），信息熵是 2 bit （当选择的参照物有两个等概率事件）：

p_A = \frac{1}{4}\\
p_B = \frac{1}{4}\\
p_C = \frac{1}{4}\\
p_D = \frac{1}{4}\\
h = log_2(4) = 2 \mathrm{bit}\\

我们也可以说，用两个硬币可以猜出选择题答案。X 事件发生的信息熵：

\mathrm{H}(X) = log_2(\text{等概率事件个数}) = \text{相对于两个等概率事件的信息熵}

如果告诉你正确答案有 1/2 概率是 A，则概率分别为：

p_A = \frac{1}{2}\\
p_B = \frac{1}{6}\\
p_C = \frac{1}{6}\\
p_D = \frac{1}{6}\\

那么我们怎么计算信息熵，自然的我们可以用概率期望的方法

\mathrm{H}(X) = p_A(log_2(?)) + p_B(log_2(?)) + p_C(log_2(?)) + p_D(log_2(?))

那等概率事件数量怎么算？概率的倒数就是等概率事件的个数：

\begin{aligned}
\mathrm{H}(\text{正确答案有 1/2 概率是 A}) &=p_A(log_2(1/p_A)) + p_B(log_2(1/p_B)) + p_C(log_2(1/p_C)) + p_D(log_2(1/p_D)) \\
&= \frac{1}{2} \times 1 + \frac{1}{6} \times (2.58) + \frac{1}{6} \times (2.58)+ \frac{1}{6} \times (2.58) \\
&= 1.79
\end{aligned}

因此我们可以得到 X 事件发生的信息熵的公式：

\mathrm{H}(X) = \sum_i^np_ilog_2(1/p_i) = -\sum_i^np_ilog_2(p_i)

那么「正确答案是 A 」信息熵是？

\begin{aligned}
\mathrm{H}(\text{正确答案是 A}) &= -\sum_i^np_ilog_2(p_i) \\
&= -1 \times 0 + 0 + 0 +0 \\
&= 0
\end{aligned}

同时，我们知道由「正确答案在ABCD之中」到被告知「正确答案有 1/2 概率是 A」的信息量是：

2 - 1.79 = 0.21 \mathrm{bit}

那么「正确答案是 A 」的信息量是？

获取信息就是消除信息熵（不确定性）

信息量 = 信息熵（不确定性）减少量

https://zhuanlan.zhihu.com/p/26486223

条件熵

熵：表示随机变量的不确定性。

条件熵：在一个条件下，随机变量的不确定性。

信息增益：熵 – 条件熵。表示在一个条件下，信息不确定性减少的程度。

在概率论中有条件概率的概念，将条件概率扩展到信息论中，就可以得到条件熵。如果两个随机变量之间具有相关性，那么在已知其中一个随机变量 X 的情况下，另一个随机变量 Y 的信息熵就会减小。

条件熵 H(Y|X) 表示的是在已知随机变量 X 的条件下另一个随机变量 Y 的不确定性，也就是在给定 X 时，根据 Y 的条件概率计算出的熵再对 X 求解数学期望：

\begin{aligned}

H(Y \mid X) 
&=\sum_{i=1}^{n} p\left(x_{i}\right) H\left(Y \mid X=x_{i}\right) \\
&=-\sum_{i=1}^{n} p\left(x_{i}\right) \sum_{j=1}^{m} p\left(y_{j} \mid x_{i}\right) \log _{2} p\left(y_{j} \mid x_{i}\right) \\
&=-\sum_{i=1}^{n} \sum_{j=1}^{m} p\left(x_{i}, y_{j}\right) \log _{2} p\left(y_{j} \mid x_{i}\right)

\end{aligned}

相对于按照 X 的不同取值进行分类，比如我们把书籍按文科和理工科分成两个书架，尽管这两堆内部依然是随机的，我们还是能根据这个信息更好地判断书籍的分布位置，但这也简化了我们的问题，使得不确定性下降。

附近学校好，买的人有5个，不买的个数为6个；

附近学校不好，买的人有1个，不买的个数为8个；

条件熵为：

\begin{aligned}
&H(Y|X = \text{附近学校好}) = -\frac{5}{11}log_2(\frac{5}{11}) -\frac{6}{11}log_2(\frac{6}{11}) = 0.99 
\\
&H(Y|X=\text{附近学校不好}) =  -\frac{1}{9}log_2(\frac{1}{9}) -\frac{8}{9}log_2(\frac{8}{9}) = 0.5\\
\\
&P(X = \text{附近学校好}) =   \frac{11}{20}\\
&P(X =\text{附近学校不好}) = \frac{9}{20}
\\ \\
&H(Y|X) =   \frac{11}{20} \times  0.99  +  \frac{9}{20} \times 0.5 = 0.77
\end{aligned}

买的人有 6 个，不买的 14 个。于是我们知道信息熵为：

H(Y) = -\frac{6}{20}log_2(\frac{6}{20}) -\frac{14}{20}log_2(\frac{14}{20}) = 0.88

信息增益 = 信息熵 – 条件熵

H(D) = H(Y) - H(Y|X)

于是得出信息增益是：

H(D_1) = 0.88 - 0.77 = 0.11

Dating-with-Machine-Learning

这个春天一起机器学习 —— Python 数据分析课。

是 Dating-with-python-this-winter 的后续。

机器学习与数据分析

目标

计划总共 14 周，每周两节理论课，

学会机器学习算法基本思想，为模型选择恰当的假设，编写机器学习代码。

前置条件：学过概率统计、线性代数。

课程表

课程文档可以下载文件，用 Typora 在电脑打开。或者上知乎 机器漫游指南。

All Cheat Sheets 是几个著名数据分析框架的“作弊表”，有些代码不会写的可以看看就上手。

课程结构

1. 机器学习入门
2. 模型评价与选择：
   1. 经验误差与结构误差
   2. 评价指标
   3. 比较检验与假设检验
3. 学习理论：频率派和贝叶斯派
4. 学习理论：PAC 和 VC维
5. 数据预处理方法
6. 学习模型：
   1. 最小二乘法线性回归——线性模型
   2. KNN 算法
   3. 决策树
   4. 支持向量机
   5. 梯度提升与随机森林
   6. 贝叶斯分类器
   7. 马尔可夫随机场

参考书目

"Learning from Data"

《统计学习方法》李航

《机器学习》周志华（西瓜书）以及《如何使用本书》

《南瓜书PumpkinBook》 辅助西瓜书

《机器学习实战：基于Scikit-Learn、Keras和 TensorFlow（原书第2版）》

为什么线性模型是可行的？

几何意义——生成子空间

假设有 N 个实例 \mathbf{x_i} = (x_0, x_1, x_2, \cdots, x_n) ， x_i(i \neq 0) 代表了实例的第 i 个属性， x_0 = 1 。

那么下列矩阵就变成了 N×(n+1) 维的数据矩阵 \mathbf{X} ，它的每一行表示同一个样本的不同属性，每一列表示不同样本中的相同属性。

\begin{aligned}
\left [  \begin{aligned} --\mathbf{x_1^\top} -- \\ --\mathbf{x_2^\top} --\\ --\dots -- \\ --\mathbf{x_N^\top}-- \\ \end{aligned} \right ]  \end{aligned}

模型是 h(\mathbf{x})=\mathbf{w}^\top \mathbf{x}=\sum_{i=0}^{n} w_{i} \cdot x_{i} ，最优解为 \mathbf{w^*} = (\mathbf{X^\top} \mathbf{X} )^{-1} \mathbf{X^\top} \mathbf{y} = X^\dagger \mathbf{y} ，代入 E_\mathrm{in}(\mathbf{w}) = \frac{1}{N} \| \mathbf{X}\mathbf{w} - \mathbf{y} \| ^2 ，得：

\begin{aligned}
E_\mathrm{in}(\mathbf{w^*}) & = \frac{1}{N}  \| \mathbf{X}\mathbf{X}^\dagger \mathbf{y}  - \mathbf{y} \| ^2 \\
& = \frac{1}{N}  \|    \mathbf{y} -  \mathbf{X}\mathbf{X}^\dagger \mathbf{y} \| ^2 \\
& = \frac{1}{N}  \|    (\mathbf{I} - \mathbf{X}\mathbf{X}^\dagger)  \mathbf{y} \| ^2 \\
& = \frac{1}{N}  \|    (\mathbf{I} - \mathbf{H} )  \mathbf{y} \| ^2
\end{aligned}

注：投影矩阵 P ∈ \mathbb{R}_{n×n} 是正交投影矩阵的充要条件 P^\top = P

如果待拟合数据任意两个属性都线性无关的话， \mathbf{X} 就可以看成一个由它的所有列向量所张成的空间。

属性的数目 n 会远远小于数据的数目 N ，因此 \mathbf{X} 张成的是 N 维空间之内的 n 维生成子空间，或者叫 n 维超平面。这个超平面的每一个维度都对应着数据集的一个列向量。理想条件下，输出 \mathbf{y} 作为属性的线性组合，也应该出现在由数据属性构成的超平面上。但受噪声的影响，真正的 \mathbf{y} 是超平面之外的一个点，这时就要退而求其次，在超平面上找到离 \mathbf{y} 最近的点作为最佳的近似。

• 黄色区域表示由所有属性张成的超平面；
• 黑色向量 \mathbf{x_1} 和天蓝色向量 \mathbf{x_2} 表示输入向量；
• 红色实线 \mathbf{y} 表示真实输出，水平的红色虚线 \mathbf{\hat{y}} 表示数据的最优估计值（属性的线性组合）；
• 垂直的红色虚线表示 \mathbf{y} - \mathbf{\hat{y}} （残差），它与超平面正交。

如果我们假设 \mathbf{w^*} 是上帝所知道的规律，那 \mathbf{Xw^*} 也仍然在 \mathbf{X} 的张成空间里。但是由于噪声 \mathbf{z} 的影响使得 \mathbf{y} 出现了偏差。

\mathbf{Xw^*} + \mathbf{z} = \mathbf{y}

这就使得

\begin{aligned}
E_\mathrm{in}(\mathbf{w^*})  
&= \frac{1}{N}\sum_{i=1}^{N}(\mathbf{y} - \mathbf{\hat{y}})^2 \\
&= \frac{1}{N}\sum_{i=1}^{N}((\mathbf{I} - \mathbf{H})\mathbf{y})^2 \\
&= \frac{1}{n}\sum_{i=1}^{N}((\mathbf{I} - \mathbf{H})\mathbf{z})^2 \\
\end{aligned}

还可以看成，这里把 \mathbf{x} 变成了它的转置（虽然输出的结果没有不同）。 \mathbf{w}^\top \mathbf{x} 背后的寓意是每个包含若干输入属性和一个输出结果的样本都被视为一个整体，误差分散在不同的样本点上；而当输出被写成 \mathbf{x}^\top \mathbf{w} 时，每个单独属性在所有样本点上的取值被视为一个整体，误差分散在每个不同的属性上。

注意我们之前令 \mathbf{x} = (x_0, x_1, x_2, \cdots, x_n) ， x_i(i \neq 0) 代表了实例的第 i 个属性， x_0 = 1 ：

h(\mathbf{x})=1 \cdot w_{0}+\sum_{j=1}^{n} x_{j} \cdot w_{j}=\mathbf{x}^{\top} \mathbf{w}

概率视角——最大似然估计 MLE

高斯噪声是最复杂的噪声，我们一般认为噪声服从正态分布：

\epsilon \sim N(\mu, \sigma^2)

真实的 f() 受到噪声的影响才有了 {y} ：

y = f(\mathbf{x}) + \epsilon = \mathbf{w^\top}\mathbf{x} +  \epsilon

y 在条件下满足概率分布：

y|_{\mathbf{x_i};\mathbf{w}}  \sim  N(\mathbf{\mu} + \mathbf{w^\top}\mathbf{x}, \sigma^2)

其概率密度函数为：

f_y(y) = \frac{1}{\sqrt{2\pi} \sigma}exp({-\frac{[y - (\mathbf{\mu} + \mathbf{w^\top}\mathbf{x})]^2} {2 \sigma^2}})

根据最大似然估计得出似然函数：

L(\mathbf{w}) = \prod_{i=1}^{N}f_y(y) \\
\begin{aligned}
ln(L(\mathbf{w})) &= \sum_{i=1}^{N} \frac{1}{\sqrt{2\pi} \sigma}exp({-\frac{[y - (\mathbf{\mu} + \mathbf{w^\top}\mathbf{\tilde{x}})]^2} {2 \sigma^2}})  \\
&= \sum_{i=1}^{N} ln(\frac{1}{\sqrt{2\pi} \sigma}) - {\sum_{i=1}^{N} \frac{[y - (\mathbf{\mu} + \mathbf{w^\top}\mathbf{\tilde{x}})]^2} {2 \sigma^2}} \\
&= \sum_{i=1}^{N} ln(\frac{1}{\sqrt{2\pi} \sigma}) - {\sum_{i=1}^{N} \frac{[y - (\mathbf{\mu} + \mathbf{w^\top}\mathbf{\tilde{x}})]^2} {2 \sigma^2}}\\
&= \sum_{i=1}^{N} ln(\frac{1}{\sqrt{2\pi} \sigma}) - {\sum_{i=1}^{N} \frac{[y - ( \mathbf{w^{*\top}}\mathbf{x})]^2} {2 \sigma^2}}

\end{aligned}

最大化似然函数：

\begin{aligned}
\mathbf{w^*} &= \arg \max _{\mathbf{w^*}}L(\mathbf{w})  \\
&=  \arg \max _{\mathbf{w^*}}\sum_{i=1}^{N} ln(\frac{1}{\sqrt{2\pi} \sigma}) - {\sum_{i=1}^{N} \frac{[y - ( \mathbf{w^{*\top}}\mathbf{{x}})]^2} {2 \sigma^2}} \\
&= \arg \min _{\mathbf{w^*}} {\sum_{i=1}^{N} \frac{[y - ( \mathbf{w^{*\top}}\mathbf{x})]^2} {2 \sigma^2}} \\
&= \arg \min _{\mathbf{w^*}} {\sum_{i=1}^{N} [y - ( \mathbf{w^{*\top}}\mathbf{x})]^2 }
\end{aligned}

对数几率回归

回归就是通过输入的属性值得到一个预测值，是否可以通过一个联系函数，将预测值转化为离散值从而进行分类呢？线性几率回归正是研究这样的问题。

为了解决一个最简单的二类分类问题, 我们为每一个点定义一个值域 [0, 1] 的函数, 表示这个点分在A类或者B类中的可能性, 如果非常可能是A类, 那可能性就逼近 1 , 如果非常可能是B类, 那可能性就逼近0（相对A的可能性）。

我们引入一个对数几率函数（logistic function ，logit 函数，也叫 sigmoid 函数）来实现实数集到 [0, 1] 的映射。将预测值投影到0-1之间，从而将线性回归问题转化为二分类问题。

y = \frac{1}{1-e^z}

一个事件发生的几率(odds)是指该事件发生的概率与该事件不发生的概率的比值。如果事件发生的概率是 p ，那么该事件的几率为 \frac{p}{1-p} ，该事件的对数几率是：

logit(y) = \ln \frac{y}{1-y}

由我们输出的得到：

\ln \frac{p(y=1 \mid \mathbf{x})}{p(y=0 \mid \mathbf{x})}=\mathbf{w}^{\mathrm{T}} \mathbf{x}+b

\begin{array}{l}
p(y=1 \mid \mathbf{x})=\frac{e^{\mathbf{w}^{\mathrm{T}} \mathbf{x}+b}}{1+e^{\mathbf{w}^{\mathrm{T}} \mathbf{x}+b}} \\
p(y=0 \mid \mathbf{x})=\frac{1}{1+e^{\mathbf{w}^{\mathrm{T}} \mathbf{x}+b}}
\end{array}

同样用 MLE：

\ell(\mathbf{w})=\sum_{i=1}^{m} \ln  p\left(y_{i} \mid \mathbf{x}_{i} ; \mathbf{w}\right)

会得到所谓的误差函数，也叫做交叉熵：

E_\mathrm{in} = \ln(1 + \exp(-y\mathbf{w^\top x}))

广义线性模型

考虑所有 y 的衍生物的情形，就得到了“广义的线性模型”（generalized linear model），其中， g 称为联系函数（link function）。

y=g^{-1}\left(\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x}+b\right)

之前的对数几率回归就是代入了 g(c) = \ln \frac{c}{1-c}

线性回归

回归问题

之前我们说，当预测值为连续值时，称为“回归问题”，离散值时为“分类问题”。但回归是什么，为什么要叫回归？

线性回归问题就是试图学到一个线性模型尽可能准确地预测新样本的输出值。

输入单一属性的问题，比如通过年龄数据预测一个人身高，输入的属性只有一个，即年龄，往往我们会先得到一系列的有标记数据，例如：[15岁，170cm] …… [20岁，175cm]。

输入多属性的问题，比如预测一个人的收入，输入的属性值就不止一个了，例如：（学历，年龄，性别，颜值，身高，体重）–>15k。回归问题就是要根据这些属性，预测新样本中人的收入。

回归的来源：生物统计学家高尔顿研究父母身高和子女身高时的发现。父亲身高和儿子身高呈正相关关系。而在正相关关系背后还有另一个现象：矮个子父亲的儿子更可能比父亲高；而高个子父亲的儿子更可能比父亲矮。高尔顿对此研究后得出的解释是自然界有一种约束力，使人类身高在一定时期是相对稳定的。如果父母身高（或矮了），其子女比他们更高（矮），则人类身材将向高、矮两个极端分化。自然界不这样做，它让身高有一种回归到中心的作用。

他当时给出了一个回归的式子，y 和 x 分别代表以英寸为单位的子代和父代的身高：

y = 3.78+0.516 x

即使父母的身高都很高，其子女不见得会比父母高，而可能会衰退（regression）（回归）至平均身高的倾向。虽然之后的x 与 y 变量之间并不总是具有“衰退”（回归）关系，但是为了纪念高尔顿这位伟大的统计学家，“线性回归”这一名称就保留了下来。

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 20, 20)
y = 3.78+ 0.516 * x
trainning_set = y + np.random.randn(y.shape[-1]) * 2
plt.plot(x, trainning_set, 'gx')

coeff = np.polyfit(x, trainning_set, 1)
poly1 = np.polyval(coeff, x)
plt.plot(x, poly1, 'b')

plt.show()

线性模型假设

假设内容：输入变量 x 的各个属性（分量）在对应的输出变量 y 中有不同的权值，或者说，输入变量的各个分量的线性组合来拟合 y 值。

\mathbf{x} = (x_1, x_2, \cdots, x_n) 是一个实例， x_i 代表了实例在第 i 个属性上的取值。我们通常令 x_0 = 1 ：

h(\mathbf{x})=\mathbf{w}^{T} \mathbf{x}=\sum_{i=0}^{n} w_{i} \cdot x_{i}

线性模型最大优点不是计算方便，而是易于解释。一些 SOTA(state of the art) 的模型里面也经常组合使用线性模型。

我们确定了模型假设，那么接下来就是确定模型的参数。

在训练集上确定系数 w_i 时，预测输出 h(x) 和真实输出 y 之间的误差是关注的核心指标。

在线性回归中，我们常常以均方误差 (MSE) 来作为模型误差。当线性回归的模型为二维平面上的直线时，均方误差就是预测输出和真实输出之间的欧几里得距离，也就是向量长度（ 或者说向量的 L2 范数）。而以使均方误差取得最小值为目标的模型求解方法就是最小二乘法。平方则是为了得到证书，因此它可以刻画样本点与直线之间的距离。

所以线性模型的泛化误差 E_{\mathrm{out}}(h) ，其中 (\mathbf{x}, y) 是未知的样本：

E_{\mathrm{out}}(h)=\mathbb{E}\left[(h(\mathbf{x}) - y)^{2}\right]

经验误差就是：

E_{\mathrm{in}}(h)=\frac{1}{N} \sum_{i=1}^{N}\left(h\left(\mathbf{x}_{i}\right) - y_{i}\right)^{2}

于是我们得到了最终目标 \mathbf{w}^* ：

\begin{array}{l}
\mathbf{w}^{*}=  \underset{\mathbf{w}}{\arg \min } E_{\mathrm{in}}(h)
\end{array}

式中每个 x_i 代表训练集中的一个样本。

求偏导以得出最值，粗体是向量或矩阵：

\begin{aligned}

E_\mathrm{in}(h) 
& = E_\mathrm{in}(\mathbf{w}) \\
& = \frac{1}{N}\sum^N_{i=1}(\mathbf{w^\top} \mathbf{x_i}  - {y_i})^2 \\
&=  \frac{1}{N}\sum^N_{i=1}(\mathbf{x_i^\top} \mathbf{w}   - {y_i})^2 \\
&=   \frac{1}{N} \left \|
\begin{aligned}
&\mathbf{x_1^\top} \mathbf{w} -   {y_1} \\
&\mathbf{x_2^\top} \mathbf{w} -   {y_2} \\
&\cdots \\
&\mathbf{x_N^\top} \mathbf{w} -   {y_N} \\
\end{aligned}
\right \| ^2 \\

&=   \frac{1}{N} \left \|
\begin{aligned}
\left [  \begin{aligned} \mathbf{x_1^\top} \\ \mathbf{x_2^\top} \\ \cdots \\ \mathbf{x_N^\top} \\ \end{aligned} \right ]  \mathbf{w} - 
\left [ \begin{aligned}  {y_1} \\ {y_2} \\ \cdots \\ {y_N} \\ \end{aligned} \right ] \end{aligned}
\right \| ^2 \\ 
& = \frac{1}{N}  \| \mathbf{X}\mathbf{w} - \mathbf{y} \| ^2

\end{aligned}

目标变成了：

\begin{aligned} 
\mathbf{w}^{*} &=\underset{\mathbf{w}}{\arg \min }  \frac{1}{N}  \| \mathbf{X}\mathbf{w} - \mathbf{y} \| ^2 \\
\\ &=\underset{\mathbf{w}}{\arg \min }   \frac{1}{N} \left(  \mathbf{w^\top} \mathbf{X^\top} \mathbf{X} \mathbf{w} - 2 \mathbf{w^\top} \mathbf{X^\top} \mathbf{y}  + \mathbf{y^\top}\mathbf{y}
\right)
\end{aligned}

求偏导：

\frac{\partial E_\mathrm{in}(\mathbf{w})} {\partial \mathbf{w}} = 
2 \mathbf{X^\top} \mathbf{X} \mathbf{w} - 2\mathbf{X^\top} \mathbf{y}

令其为0 ，考虑到矩阵不可逆（伪逆），得出：

\mathbf{w^*}  =  (\mathbf{X^\top} \mathbf{X} )^{-1} \mathbf{X^\top} \mathbf{y}
= X^\dagger \mathbf{y}

在单变量线性回归任务中，最小二乘法的作用就是找到一条直线，使所有样本到直线的欧式距离之和最小。说到这里，问题就来了：凭什么使均方误差最小化的参数就是和训练样本匹配的最优模型呢？