统计机器学习03：线性模型
本页统计信息
前面说到机器学习中的任务主要是分类和回归，分类如前面的贝叶斯定理所展现的那样是输出一个类别(有限的离散数值)，而回归的任务则是输出一个real value
线性回归的基本思路
对于n维的样本\(\boldsymbol x=[x_1,x_2,…,x_n]\)找到一系列线性函数\(\omega=[\omega_1,\omega_2,…,\omega_n]\)使得\(f(x)=\omega^T\boldsymbol x+b\)，这个问题可以变形为n+1维使得\(f(x)=\omega^T\boldsymbol x\)其中\(\boldsymbol x=[1,x_1,x_2,…,x_n]\)而\(\omega=[b,\omega_1,\omega_2,…,\omega_n]\)，这样一来表达式的形式更加统一了，而我们的目标就是训练出这个函数f，使得对于训练数据\((x_i,y_i)\)，有\(f(x_i)=y_i\)，我们要求解的对象就是这个n+1维的向量\(\omega\)
线性回归的loss函数和解
线性回归常见的loss函数定义是最小平方误差(MSE)
\[
MSE=\frac{1}{n}\sum\limits_{i=1}\limits^{n}(y_i-f(x_i,\omega))^2
\]
也会使用残差平方和(RSS)，其形式如下：
\[
J_n(\alpha)=\sum\limits_{i=1}\limits^{n}(y_i-f(x_i,\omega))^2=(\boldsymbol y-X^T\omega)^T(\boldsymbol y-X^T\omega)
\]
我们对RSS的表达式求梯度有：
\[
\nabla J_n(\alpha)=-2X(\boldsymbol y-X^T\alpha)=0
\]
因此可以得到使得RSS最小的线性回归的解是:
\[
\omega=(XX^T)^{-1}Xy
\]
我们需要注意到，每一个样本x是\(d\times 1\)维的向量，因此X是\(d\times n\)维的矩阵，而y也是\(n\times 1\)维的向量，所以当样本数n小于特征数d的时候，\(XX^T\)是不满秩的，求不出逆矩阵，此时的线性回归有多个解
其实这也很好理解，因为要求解的是n+1维的向量\(\omega\)，有n+1个变量，因此至少需要n+1个样本才能确定n+1个参数，出现上述情况的时候所有可能的解都可以使得均方误差最小化，此时可以考虑引入正则化项来筛选出需要的结果。
线性回归的统计模型
真实情况下的数据样本往往会有噪声，比如
\[
y=f(\boldsymbol x,\omega)+\epsilon
\]
其中\(\epsilon\)是一个随机噪声，服从\(N(0,\sigma^2)\)的正态分布,此时可以通过极大似然法来估计\(\omega\)，定义
\[
P(y|\boldsymbol x,\omega,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}\exp[-\frac{1}{2\sigma^2}(y-f(\boldsymbol x,\omega))^2]
\]
根据极大似然法有
\[
L(D,\omega,\sigma)=\prod_{i=1}\limits^nP(y_i|\boldsymbol x_i,\omega,\sigma)
\]
我们的求解目标变成了：
\[
\omega=\arg\max L(D,\omega,\sigma)= \arg \max\prod_{i=1}\limits^nP(y_i|\boldsymbol x_i,\omega,\sigma)
\]
取对数似然之后有
\[
l(D,\omega,\sigma)=-\frac{1}{2\sigma^2}\sum_{i=1}\limits^{n}(y_i-f(x_i,\omega)^2)+c(\sigma)
\]
到这一步为止我们又回到了RSS，因此解的表达式依然和上面推出的是一样的。
岭回归 Ridge Regression
我们发现普通的线性回归很容易出现overfitting的情况，比如一些系数\(w_i\)的值非常极端，仿佛就是为了拟合数据集中的点而生的，对测试集中的数据表现就非常差。为了控制系数的size，让表达式看起来更像是阳间的求解结果，我们可以引入正则化的方法。
\[
\omega^{*}=\arg\min \sum_{i=1}\limits^{n}(y_i-w^Tx_i)^2+\lambda\sum_{j=1}\limits^{d}w_j^2
\]
这实际上就是拉格朗日算子，则我们跟之前一样可以将其写成矩阵形式：
\[
(\boldsymbol y-X^T\omega)^T(\boldsymbol y-X^T\omega)+\lambda\omega^T\omega
\]
对其求梯度可以得到
\[
\nabla J_n(\alpha)=-2X(\boldsymbol y-X^T\alpha)+2\lambda\omega=0
\]
则其最终的解就是：
\[
\omega^{*}=(XX^T+\lambda I)^{-1}Xy
\]
这就是岭回归(Ridge Regression)的方法，其中参数\(\lambda\)是可以自己确定的，我们可以保证矩阵\(XX^T+\lambda I\)是满秩的矩阵。
贝叶斯线性回归
现在我们重新考虑实际样本中可能会出现的噪声，即\(y=f(\boldsymbol x,\omega)+\epsilon\)，其中\(\epsilon\)服从\(N(0,\sigma^2)\)的正态分布。
根据贝叶斯定理可以得到：
\[P(\omega|y,x,\sigma)=\frac{P(y|\omega,x,\sigma)P(\omega|x,\sigma)}{P(y|x,\sigma)}
\]
即posterior正比于prior和likelihood的乘积，即\(\ln(posterior)\propto \ln(likelihood)\times\ln(prior)\)，而在线性回归问题中，我们已经
知道了likelihood就是
\[l(D,\omega,\sigma)=-\frac{1}{2\sigma^2}\sum_{i=1}\limits^{n}(y_i-f(x_i,
\omega)^2)+c(\sigma)
\]
我们可以用如下方法选择：
\[
\begin{aligned}
p(\omega) =N(w|0,\lambda^{-1}I)&=\frac{1}{(2\pi)^{\frac{d}{2}}|\lambda^{-1}I|^{\frac{1}{2}}}\exp^{-\frac{1}{2}\omega^T(\lambda^{-1}I)\omega}\\
\ln(p(\omega)) &=-\frac{\lambda}{2}\omega^T\omega+c\end{aligned}
\]
因此在贝叶斯理论下的岭回归模型需要优化的目标可以等价于：
\[
-\frac{1}{2\sigma^2}\sum_{i=1}\limits^{n}(y_i-f(x_i,\omega)^2)+c(\sigma)-\frac{\lambda}{2}\omega^T\omega+c
\]
逻辑回归 Logistic Regression
基本概念
逻辑回归往往选择通过一个sigmod函数将样本映射到某个区间上，以此来估计样本属于某种类别的概率从而达到分类的目的。我们经常选用的sigmod函数是：
\[
y=\sigma(z)=\frac{1}{1+e^{-z}}
\]
比如对于二分类问题，可以将样本映射到区间(-1,1)上，此时如果计算结果为正数则说明样本属于正例，反之就是反例，可以表示为：
\[
\begin{aligned}
P(y_i=1|x_i,\omega)=\sigma(\omega^Tx_i)=\frac{1}{1+e^{-\omega^Tx_i}}\\
P(y_i=-1|x_i,\omega)=1-\sigma(\omega^Tx_i)=\frac{1}{1+e^{\omega^Tx_i}}\end{aligned}
\]
上面的两个式子也可以统一写为：
\[
\begin{aligned}
P(y_i=1|x_i,\omega)=\sigma(y_i\omega^Tx_i)=\frac{1}{1+e^{-y_i\omega^Tx_i}}\end{aligned}
\]
参数估计
我们可以用极大似然法来估计参数\(\omega\)，依然用D表示数据集，具体的过程如下所示：
\[
\begin{aligned}
P(D) & =\prod_{i\in I}\sigma(y_i\omega^Tx_i)\\
l(P(D)) & =\sum_{i\in I}\ln(\sigma(y_i\omega^Tx_i))=-\sum_{i\in I}\ln(1+e^{y_i\omega^Tx_i})
\end{aligned}
\]
因此我们可以将逻辑回归的极大似然法参数估计的loss函数定义成：
\[
E(\omega)=\sum_{i\in I} \ln(1+e^{-y_i\omega^Tx_i})
\]
对于一个二分类问题，我们如果用0和1而不是+1和-1来代表两种分类，那么上面的表达式又可以写为：
\[
\begin{aligned}
E(\omega)&=\sum_{i\in I\cap y_i=1}\ln(1+e^{-\omega^Tx_i}) + \sum_{i\in I\cap y_i=0}\ln(1+e^{\omega^Tx_i})\\
&=\sum_{i\in I\cap y_i=1}\ln (e^{\omega^Tx_i})(1+e^{\omega^Tx_i}) + \sum_{i\in I\cap y_i=0}\ln (1+e^{-\omega^Tx_i})\\
&=\sum_{i\in I} \ln(1+e^{\omega^Tx_i})-\sum_{i\in I \cap y_i=1}e^{\omega^Tx_i} \\
&=\sum_{i\in I} (-y_i\omega^Tx_i+\ln(1+e^{\omega^Tx_i}))
\end{aligned}
\]
我们可以证明\(E(\omega)\)是一个关于\(w\)的凸函数，根据凸函数的可加性，我们只需要证明\(-y_i\omega^Tx_i+\ln(1+e^{\omega^Tx_i})\)是关于\(w\)的凸函数，我们令\(g(\omega)=-y_i\omega^Tx_i+\ln(1+e^{\omega^Tx_i})\)则对其求一阶梯度可以得到：
\[
\frac{\partial g(\omega)}{\partial\omega} = -y_ix_i+\frac{x_ie^{\omega^Tx_i}}{1+e^{\omega^Tx_i}}
\]
能得到这个结果是因为我们有这样一个结论：对于一个n维的向量\(\omega\),我们有\(\frac{\partial\omega^T}{\partial\omega}=I_{n}\)
进一步地，我们对上面求的一阶梯度再求二阶梯度，可以得到：
\[
\frac{\partial^2g(\omega)}{\partial\omega^2}=\frac{x_i^2e^{w^Tx_i}}{(1+e^{w^Tx_i})^2}\geq 0
\]
因此我们证明了损失函数是一个凸函数，因此可以用梯度下降的方法求得其最优解，即：
\[
\omega^{*}=\arg\min_{\omega} E(\omega)
\]
根据上面求得的一阶梯度，可以得到基于梯度下降法的逻辑回归参数求解迭代方程：
\[
\begin{aligned}
\omega_{i+1}&=\omega_{i}-\eta(i)\sum_{i\in I}(-y_ix_i+\frac{x_ie^{\omega_i^Tx_i}}{1+e^{\omega_i^Tx_i}})\\
&=\omega_i+\eta(i)\sum x_i(\frac{1}{1+e^{-\omega_i^Tx_i}}-y_i)\\
&=\omega_i+\eta(i)X(\sigma(\omega_i, X)-y)
\end{aligned}
\]
其中\(\sigma(x)\)是sigmod函数，而\(\eta(i)\)是自己选择的学习率，一半是一个比较小的数字，并且应该不断减小。