【机器学习】1. 广义线性模型

2023-03-10 14:35 由江景景景页发表于 #后端开发

1. 广义线性模型

1.1 从简单线性回归开始

机器学习的任务是从已知的数据中提取知识, 从而完成对新数据的识别与预测, 即应用知识. 但是, 数据本身不会给予研究者想要的答案, 大多数时候人们很难通过观察或者简单的计算提取有用的结论. 为了从历史数据中整合知识, 人们提出了"模型"这一概念, 认为数据是从模型中生成的, 遵循一定的规律, 如果规律是可知的, 那么机器从数据中学习规律就是可能的.

由于数据不会说话, 最初的模型由人们假设(assumption)而来. 下面用大家熟知的简单线性回归为例, 描述一个数据生成的过程. 假设数据为\(\{(y_i,x_i)\}_{i=1}^{n}\), 其中\(x_i\)是一个向量, 描述样本点\(i\)的特征(features). 简单线性回归假设所有的样本点由下面的模型生成:

\[y=\theta^Tx+\epsilon \]

这里\(\theta\)是一个与\(x\)同型的向量, 代表了模型的可知部分, \(\epsilon\)是一个不可知的正态零均值随机误差, 方差为\(\sigma^2\). 上述模型也可以写成下面的概率分布形式:

\[y\sim N(\theta^Tx,\sigma^2). \]

我们获得了从这个模型中生成的数据集\(\{(y_i,x_i)\}_{i=1}^{n}\), 希望从中学习到关于可知部分\(\theta\)的知识, 因为知道了\(\theta\), 就可以找到数据生成的过程, 从而对任何新观测的特征\(x_{n+1}\), 即使不能知道对应的随机误差, 也能够大概了解标签\(y_{n+1}\)的期望\(\mathrm{E}[y_{n+1}]\), 这就是把从历史数据中学到的知识应用于新数据的场景.

现在的目标是尽可能准确地估计一个\(\theta\)的值\(\hat{\theta}\)得到估计模型, 因为真实的\(\theta\)是不会直接呈现给我们的, 只能从历史数据中估计. 对于每一个估计的\(\hat{\theta}\), 得到一个估计的概率分布\(N(\hat{\theta}^Tx,\sigma^2)\), 即我们视为数据集是从这个分布中产生的, 在期望意义下, \(h_{\hat\theta}(x)=\mathrm{E}[N(\hat{\theta}^Tx,\sigma^2)]=\hat{\theta}^Tx\)是从这个概率分布中抽样最应该出现的值.

如果模型准确, \(h_\hat{\theta}(x)\)和真实数据\(y\)的差异肯定尽可能小, 因此只需要想办法刻画模型结果和真实数据的差异并最小化之. 对两个在\(\mathbb{R}^{n}\)上取值的向量, 衡量其差异的方法有很多, 选择均方误差作为其差异的度量, 要最小化两个向量之间的均方误差, 意味着我们的目标(objective)是

\[\hat\theta=\min_{\hat\theta}(h_\hat\theta(x)-y)^T(h_{\hat\theta}(x)-y). \]

优化不属于我们讨论的范畴, 我们只需要知道梯度下降法可以优化这个目标函数, 得到一个\(\hat\theta\)即可. 这样, 我们就得到了估计模型.

但是, 简单线性回归模型显得过于简单了. 如果我们想获得一个更加有应用价值的模型, 至少有下面几点是可以考虑的:

在数据生成机制上, 没有任何的非线性算子使得数据的生成机制被我们描述得过于理想化, 但过多的非线性也会增大估计难度. 因此可以考虑保留核心的线性生成部分\(\theta^Tx\), 但对结果施加一个非线性函数.
模型可能不由正态分布所生成, 可以尝试其他分布.
我们也许感兴趣的不是随机数\(y\)本身, 而是其某种变换\(T(y)\).

考虑以上三点, 就从简单线性回归模型进入了广义线性模型.

1.2 指数族

现在先考虑前两点, 尝试某种其他的分布\(\mathrm{Distribution}(\eta)\), 其中核心未知参数是\(\theta^Tx\)的某个变换, 即

\[y\sim \mathrm{Distribution}(a(\theta^Tx)). \]

假定我们知道分布的具体形式, 也获得了数据集\((y,x)\), 那么对于任何\(\hat{\theta}\), 如果数据来源于模型\(\mathrm{Distribution}(a(\hat\theta^Tx))\), 那么此分布的期望和真实数据\(y\)必定差异较小(这里我们用"距离"来表示, 以区别于一般的作差), 也就是说我们的目标依然是

\[\min_{\hat{\theta}} \mathrm{Distance}(y,\mathrm{E}[\mathrm{Distribution}(a(\hat\theta^Tx))]). \]

在简单线性回归中, 距离函数被我们选择为均方差, 现在我们且不论距离的选择. 现在摆在面前的问题是如何计算分布的期望, 如果有可能的话, 最简单的期望形式必然是\(a(\theta^Tx)\)的某个映射, 这让我们不能对所选分布太过随意. 有一类特殊的分布满足我们的需求: 指数族(exponential family). 由于指数族定义不依赖于模型, 先用\(\eta\)代替\(\theta^Tx\).

如果某个参数分布族的密度函数满足下面的形式:

\[p(y;\eta)=b(y)\exp(\eta T(y)-a(\eta)), \]

就称此分布族是指数分布族, 这里\(a(\cdot)\), \(b(\cdot)\)和\(T(\cdot)\)都是可微函数. 我们关注它的期望:

\[\mathrm{E}[p(y;\eta)]=\int_{\mathbb{R}}yp(y;\eta)\mathrm{d}y=\int_{\mathbb{R}} yb(y)\exp(\eta T(y)-a(\eta))\mathrm{d}y, \]

直接求期望较困难, 注意到对分布族中的任意概率分布, 密度函数的积分为\(1\), 因此我们可以尝试将指数族密度函数的积分对参数\(\eta\)求导, 其导数值应恒为\(0\). 指数族积分和求导可交换, 为计算提供了便利([2]中给出了不使用可交换性条件下的证明):

\[\begin{aligned} \frac{\partial }{\partial \eta}\int_{\mathbb{R}} p(y;\eta)\mathrm{d}y&=\int_{\mathbb{R}}\frac{\partial b(y)\exp(\eta T(y)-a(\eta))}{\partial \eta}\mathrm{d}y\\ &=\int_{\mathbb{R}}(T(y)-a'(\eta))b(y)\exp(\eta T(y)-a(\eta))\mathrm{d}y\\ &=T(\mathrm{E}[p(y;\eta)])-a'(\eta)\\ &=0, \end{aligned} \]

这样我们求出了\(\mathrm{E}[p(y;\eta)]\)的一个变换的值, 如果\(T(\cdot)\)是可逆函数, 那么\(\mathrm{E}[p(y;\eta)]=T^{-1}(a'(\eta))\); 特别地如果\(T(\cdot)\)是恒等变换, 那么\(\mathrm{E}[p(y;\eta)]=a'(\eta)\). 这里我们得到了指数族的期望, 证明指数族正是我们所关注的分布族.

现在回到机器学习任务上, 假如我们选择的分布是指数族分布\(\mathrm{ExponentialFamily}(\theta^Tx)\), 并取\(T(\cdot)\)为恒等映射, 那么根据在简单线性回归一节中所讨论的, 我们只需要:

对每一个估计的\({\theta}\), 代入指数族中, 得到期望观测向量\(a'(\theta^Tx)\).
选择一个合适的距离函数, 最小化\(\mathrm{Distance}(a'(\theta^Tx),y)\).
如果最小化由迭代完成, 就重复上述两个步骤; 否则用其他方式最小化.

1.3 广义线性回归实例

现在回顾简单线性回归模型, 由\(N(\theta^Tx,\sigma^2)\)也就是\(N(\eta,\sigma^2)\)生成, 我们考虑将正态分布化为指数族的形式:

\[p(y;\eta)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{y^2-2\eta y+\eta^2}{2\sigma^2} \right)=\frac{\exp(-\frac{y^2}{2\sigma^2})}{\sqrt{2\pi\sigma^2}}\cdot\exp\left(\frac{\eta\cdot y-\eta^2}{\sigma^2} \right), \]

这里存在另一个参数\(\sigma^2\), 这关系到指数族一个更广泛的形式:

\[p(y;\eta,\tau)=b(y,\tau)\exp\left(\frac{\eta^TT(y)-a(\eta)}{c(\tau)} \right). \]

这里\(\eta=\theta^Tx\)是向量, \(\theta\)是参数矩阵, \(T(y)\)是一个向量函数, 这种表述下文还会提到. 上式中\(\tau\)称为扩散系数, 也就是正态分布中的\(\sigma^2\). 另外, 指数分布族中包含众多分布: 正态分布, 两点分布, 二项分布, 泊松分布, Gamma分布等.

以两点分布\(B(1,p)\)为例, 它常常被视为二分类任务的生成模型. 要注意不能直接用\(\eta\)替换\(p\), 至少因为\(p\)的取值范围只有\([0,1]\)而\(\eta=\theta^Tx\in\mathbb{R}\), 具体\(\eta\)应该如何表现还应当先将两点分布的概率函数(在这里为分布列)改写为指数族的形式.

\[p(y;p)=p^{y}(1-p)^{1-y}=\exp\left(y\log p+(1-y)\log(1-p) \right)=\exp\left(y\cdot\log\frac{p}{1-p}+\log(1-p) \right). \]

这个密度形式与指数族的标准形式不同, 但只要施加变换

\[\eta=\log\frac{p}{1-p}, \]

同时有

\[p=\frac{e^{\eta}}{1+e^\eta}, \]

就可以把密度变为

\[p(y;p)=\tilde p(y;\eta)=\exp\left(y\cdot\eta-\log(1+e^{\eta}) \right), \]

变成指数族标准形式, 且\(a(\eta)=\log(1+e^{\eta})\). 对比前面给出的结论, 从\(B(1,p)\)中抽出的样本\(y\)应当与\(a'(\eta)\)一致, 即最小化

\[\mathrm{Distance}(y,a'(\eta))=\mathrm{Distance}\left(y,\frac{e^{\theta^Tx}}{1+e^{\theta^Tx}} \right)=\mathrm{Distance}\left(y,\mathrm{sigmoid}(\theta^Tx) \right). \]

这就是二分类器中, 用线性函数sigmoid变换作概率估计的由来. 同时, 对分类问题, 尽管均方差用作距离函数依然是可行的, 但梯度下降中均方差表现不好, 因此我们往往使用交叉熵函数作距离度量(原因这里不谈):

\[\mathrm{CrossEntropy}(y,\hat y)=-\sum_{i=1}^{n}y_i\log_2 \hat{y}_i. \]

最后, 回到上文我们提到的\(\theta\)为参数矩阵的情况, 此时对应的分布族是多项分布\(B(1;p)\), 其中\(p=(p_1,\cdots,p_k)'\), 每次生成一个单位向量\(y=(y_1,\cdots,y_k)\)其中有且仅有一个分量\(y_i=1\)(以概率\(p_i\)), 其余分量为\(0\), , 它常常被认为是多分类任务的生成模型.

尝试写出多项分布的概率函数, 注意, 一个生成\(k\)维向量的多项分布中, 未知参数只有\(k-1\)个, 这是因为\(\sum_ip_i=1\)且\(\sum_iy_i=1\).

\[\begin{aligned} p(y;p)&=p_1^{y_1}p_2^{y_2}\cdots p_k^{y_k} \\ &=\exp\left(\sum_{i=1}^{k-1}y_i\ln p_i+\left(1-\sum_{i=1}^{k-1}y_i \right)\ln\left(1-\sum_{i=1}^{k-1}p_i \right) \right)\\ &=\exp \left(\sum_{i=1}^{k-1}y_i\ln\frac{p_i}{1-\sum_{i=1}^{k-1}p_i}+\ln\left(1-\sum_{i=1}^{k-1}p_i \right) \right) \end{aligned} \]

至此, 令

\[\eta=\left(\ln\frac{p_1}{p_k},\cdots,\ln\frac{p_{k-1}}{p_k},0 \right), \]

可得\(p_i=p_ke^{\eta}(i=1,\cdots,k-1)\), 结合\(\sum_ip_i=1\)可得

\[p=\left(\frac{e^{\eta_1}}{1+\sum_{i=1}^{k-1}e^{\eta_i}},\cdots,\frac{e^{\eta_{k-1}}}{1+\sum_{i=1}^{k-1}e^{\eta_i}},\frac{1}{1+\sum_{i=1}^{k-1}e^{\eta_i}} \right), \]

上面两个繁琐的式子实际上就是\(\eta_i=\ln(p_{k-1}/p_k)\)以及\(p_i=e^{\eta_i}/\sum_je^{\eta_j}\). 此时密度函数为

\[p(y;p)=\tilde p(y;\eta)=\exp\left(y^T\eta-\ln \left(\frac{e^{\eta_k}}{\sum_{i=1}^{k}e^{\eta_i}} \right)\right):=\exp(y^T\eta-a(\eta)) \]

这里\(a'(\eta)\)就是对其每个分量求导, 即

\[a'(\eta)=\begin{bmatrix} \frac{e^{\eta_1}}{\sum_i e^{\eta_i}} \\ \vdots \\ \frac{e^{\eta_k}}{\sum_{i}e^{\eta_k}} \end{bmatrix}=\begin{bmatrix} \frac{\exp(\theta_1^Tx)}{\sum_i\exp(\theta_i^Tx)}\\ \vdots \\ \frac{\exp(\theta_k^Tx)}{\sum_i\exp(\theta_i^Tx)} \end{bmatrix}. \]

此时\(\theta\)是参数矩阵, \(\theta_i\)代表矩阵的第\(i\)列.

1.4 广义线性模型代码

statsmodels库中提供了单参数指数族对应的广义线性模型. 以下假设我们有一个从泊松分布中生成的模型, 这里先验证泊松分布族是指数族.

\[p(y;\lambda)=\frac{\lambda^y}{y!}e^{-\lambda}=\frac{1}{y!}\exp(y\log \lambda-\lambda)\xlongequal{\eta=\log\lambda}\frac{1}{y!}\exp(y\eta-e^{\eta}), \]

也即数据从\(P(e^{\theta^Tx})\)中生成, 生成数据的期望为\(a'(\eta)=e^{\eta}=e^{\theta^Tx}\). 现在我们先生成数据. 假设\(x\)是三维数据, \(\theta=(1, 2, 5)^T\).

import numpy as np
import statsmodels.api as sm
np.random.seed(2000)

## Generate data from poisson distribution

theta = np.array([1, 2, 5])
x = np.random.normal(size=(100, 3))
y = np.random.poisson(lam=np.exp(np.dot(x, theta))) + np.random.normal(scale=1, size=(100,))

注意, GLM()函数是不会自己添加截距项的, 这里简化处理, 我们也认为模型不含截距项. 因为是泊松分布族生成的数据, 因此要选择分布族为Poisson().

glm = sm.GLM(y, x, family = sm.families.Poisson())
res = glm.fit()
print(res.summary())

得到的结果为\(\hat{\theta}=(0.9964,1.9946,5.0010)'\), 非常接近真实值.

在GLM()函数中可以设定联结函数(link function), 每一个分布族会自带默认的联结函数, 这个函数就是连接分布参数\(\lambda\)和自然参数\(\eta=\theta^Tx\)的函数, 本例中为\(\log\), 因为\(\eta=\log\lambda\). 有关联结函数的更多信息, 可参考[3].

Reference

[1] https://cs229.stanford.edu/notes2021fall/cs229-notes1.pdf

[2] https://xg1990.com/blog/archives/304

[3] https://www.statsmodels.org/stable/glm.html

热门相关：斗神战帝网游之逆天飞扬豪门闪婚：帝少的神秘冷妻刺客之王大神你人设崩了