机器学习之广义线性模型

机器学习之广义线性模型

在机器学习中，我们常常是从线性回归和Logistics回归这两种模型入手。大多数人在学的时候时将其当成两个独立的模型去学习的，线性回归用来拟合直线，Logistics回归用来分类。实际上，这两种模型都是一个更广泛模型的特例，这就是广义线性模型（Generalized Linear Models, 简称GLM）。

• 引子：线性回归和Logistics回归
• 从概率的角度解释回归
• 指数型分布族（指数族）
• 指数型分布族的向量化写法
• 自然指数族
• 指数分散族
• 广义线性模型与指数型分布族（核心）
• 线性的体现
• 连接函数与激活函数
• 标准连接
• 广义线性模型的参数关系
• 广义线性模型的最大似然估计
• 标准连接函数下的广义线性模型最大似然估计
  • Logistics回归优化举例
  • 补充：牛顿法的简化方法之一Fisher分数法
• 回答引子的疑问，最大似然估计形势下的迭代优化
• 广义线性模型的求解（IRLS算法）
• 参考文献

引子：线性回归和Logistics回归

如果刚学完线性回归和Logistics回归，那么是否会注意到，二者的梯度更新步骤都是(虽然\(h_{\vec\theta}(\vec x^{(i)})\)的定义不同)： \[ \theta_j=\theta_j-\alpha(h_{\vec\theta}(\vec x^{(i)})-y^{(i)})x_j^{(i)}\\ h_{\vec\theta}(\vec x^{(i)})=\begin{cases} \vec{\theta}^T \vec{x},\quad线性回归\\ \frac{1}{1+e^{-\vec{\theta}^T \vec{x}}},\quad Logistics回归\end{cases} \] 其中，\(\vec\theta, \vec x^{(i)}\)分别是参数向量，第\(i\)个观测数据的向量。下标\(j\)表示第\(j\)个分量，\(\alpha\)表示更新的步长。那么他们二者为什么都有相同的更新公式呢？(只有\(h_{\vec{\theta}}(\vec{x})\)的具体表现形式不同)。第一个原因是他们本质上都可以从最大似然估计推导出来；第二个原因则是它们二者都是一种更普遍的模型的特殊情况，这个模型就是广义线性模型。

从概率的角度解释回归

如果我们详细看统计回归的过程，就能发现它有两步组成。第一步是参数估计，确定特定模型中的未知参数，即求模型的参数\(\Theta\)；第二部根据已经确定的模型与参数，预测新数据的函数值，即求\(y_{pred}\)。

我们使用参数估计的过程主要是使用MLE、MAP、Bayesian等准则（推荐文章ML, MAP, and Bayesian --- The Holy Trinity of Parameter Estimation and Data Prediction），这种方式往往是以上述某个准则推导出的概率最大的值作为参数估计的结果。因此，在第二部预测步骤时，给定参数估计结果\(\theta\)和自变量\(x\)后，通过模型得到的预测结果往往也不是实际真实值，而是真实值和某个概率分布相关的误差组合起来的结果（我们希望这个组合出来的概率分布期望就是真实值）。

回看线性回归的结果，其预测模型是一条直线，但是真实的数据点并不一定在直线上，而是以某个概率分布在预测模型周围。

以线性回归模型为例，假设回归函数为\(y=\mathbf{\theta}^T\mathbf{x}\) (\(\theta,x\)为向量), 对于每对观测结果\((x^{(i)},y^{(i)})\)，都有 \[y^{(i)}=\theta^Tx^{(i)}+\epsilon^{(i)}\] 其中 \(\epsilon\)为误差，基于一种合理的假设（中心极限定理），我们可以认为误差的分布服从正态分布(又称高斯分布)，即 \(\epsilon \sim N(0,\sigma^2)\) ，那么，我们可以认为\(y^{(i)} \sim N(\theta^Tx^{(i)},\sigma^2)\),根据正态分布的概率公式 \[P(y^{(i)}|x^{(i)};\theta)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})\]

我们可以通过最大似然估计求出模型的参数\(\mathbf{\theta}\)，进而得到线性模型\(y=\mathbf{\theta}^T\mathbf{x}\)。我们不难发现，\(\mathbf{\theta}^T\mathbf{x}\)正好是\(y\)这个高斯随机变量的期望！

线性回归这条直线\(y=f(x)\)的真实含义其实是：回归模型中，对于取特定的自变量\(x^{(i)}\)，其因变量概率的期望是\(f(x^{(i)})\)。同样，在Logistics回归中，最终的预测结果也是二项分布中取0或1的期望，其以0.5作为阈值的原因也在于此，若以大于0.5的期望取1，那么我们就认为结果是1；若以小于0.5的期望取1（取0概率大于0.5），那么我们就认为结果是0，本质上就是这么直白。

此外，也可从\(y_{pred}\)为一个统计量的角度理解。

我们之前已经说过，我们使用回归模型预测的值，其实也并不是精确值，而是预测值概率的期望。广义线性模型和概率的关系，就是我们的预测值\(y_{pred}\)的概率分布服从指数族分布。而指数族的参数通过连接函数和线性函数连接到一起，即\(g(u)=\theta^T x\)，这一点我们之后再说。下面我们先看指数族。

指数型分布族（指数族）

指数型分布族是指数分布族的推广，囊括了正态分布族、二项分布族、伽马分布族、多项分布族常见分布等等。具体定义形式如下：

指数型分布族：一个概率分布族\(\mathfrak{p}=\{p_{\theta}(x);\theta∈\varTheta\}\)可称为指数型分布族，假如\(\mathfrak{p}\)中的分布（分布列或密度函数）都可表示为如下形式： \[p_\theta(x)=h(x)c(\theta)\exp\left\{\sum_{j=1}^k c_j(\theta)T_j(x)\right\}\tag{1}\] 其中，k为自然数；\(\theta\)可以是数字，也可以是向量。分布的支撑\(\{x:p(x)>0\}\)与参数\(\theta\)无关；诸\(c(\theta),c_1(\theta),\dotsb,c_k(\theta)\)是定义在参数空间\(\varTheta\)上的函数；诸\(T_1(x),\dotsb,T_k(x)\)是\(x\)的函数，称为充分统计向量，但\(T_1(x),\dotsb,T_k(x)\)线性无关。\(h(x)\)也只是\(x\)的函数，且\(h(x)>0\)，通常是一个常数。

\(c(\theta)\)是作为归一化参数存在的，称为叫做配分函数(partition function)。 \[c(\theta)^{-1} = \int h(x) \exp\left\{\sum_{j=1}^k c_j(\theta)T_j(x)\right\} dx\] 此外，指数族还有另一种表述方式，就是将外面的\(c(\theta)\)放到指数符号中： \[p_\theta(x)=h(x)\exp\left\{\sum_{j=1}^k c_j(\theta)T_j(x)-A(\theta)\right\}\tag{2}\] 由于通常\(A(\theta)\)含有\(\log\)符号，该部分也称为“Log Partition Function”，易知\(A(\theta)=\ln c(\theta)\)。 如果我们使用向量值函数来表达指数型分布族可写为: \[p_\theta(x)=h(x)\exp\left\{\sum_{j=1}^k c_j(\theta)T_j(x)-A(\theta)\right\}\tag{3}\]

从上述定义可知，一个分布族是不是指数型分布族的关键在于其概率分布能否改写为定义中方式。

指数型分布族的向量化写法

下面我们使用向量值函数将式(4)进行进一步改造。

向量值函数：有时也称为向量函数，是一个单变量或多变量的、值域是多维向量或者无穷维向量的集合的函数。向量值函数的输入可以是一个标量或者一个向量，输出是向量，定义域的维度和值域的维度是不相关的。

对于\(\theta\)的一系列函数\(c_1(\theta),c_2(\theta),\dotsb\)和充分统计量向量\(T_1(x),T_2(x),\dotsb\)，我们写出列向量形式： \[ \mathbf{C}(\theta)=\begin{bmatrix}c_1(\theta)\\c_2(\theta)\\\vdots\\c_k(\theta)\end{bmatrix} \mathbf{T}(x)=\begin{bmatrix}T_1(x)\\T_2(x)\\\vdots\\T_k(x)\end{bmatrix} \] 那么式（3）可写成 \[ p(x;\theta)=h(x)\exp\left\{\mathbf{C}^T(\theta)\mathbf{T}(x)-A(\theta)\right\}\tag{4} \] 其中，\(\mathbf{C}(\theta),\mathbf{T}(x)\)都是向量值函数，\(h(x),A(\theta)\)都是普通函数。通常文章会把\(A(\theta)\)写成\(A(\mathbf{C}(\theta))\)的形式，这两种本质上是等价的，但是\(A(\mathbf{C}(\theta))\)的参数形式更加统一，为主流用法。由于\(\mathbf{C}(\theta)\)的计算结果本质上就是一个向量，我们可令向量值函数\(\mathbf{C(\theta)}=\eta\)，那么式（4）可表示为： \[ p(x;\eta)=h(x)\exp\left\{\eta^T\mathbf{T}(x)-A(\eta)\right\}\tag{5} \] 这就是其他资料中的常见形式。其中\(\eta=\mathbf{C}(\theta)\)，参数\(η\)通常叫做自然参数(natural parameter)或者标准参数(canonical parameter)。这里注明：\(A(\theta)\)与\(A(\eta)\)实际上是两个不同的函数，但是可以通过\(\eta=\mathbf{C}(\theta),\theta=\mathbf{C}^{-1}(\eta)\)进行互换，因此在后文对他们不做区分。此外，\(\eta,\theta\)是一一对应的，这里先不加证明地写出这个引理。

引理1：在指数族中函数\(C(\cdot)\)总是单调连续的(存在逆函数)，所以自然参数\(η\)和原始参数\(θ\)是存在一一映射关系的。 \[ \eta=\mathbf{C}(\theta)\\ \theta=\mathbf{C}^{-1}(\eta) \]

在指数型分布族中，使用标准参数\(η\)表示的公式形式称为指数族分布的标准形式(canonical form)，在标准形式下，分布的参数是\(η\)。实际上，从原始分布向指数型分布转换的过程就是将\(\theta\)转换为\(\eta\)的过程。

广义线性模型的预测值概率分布都属于指数型分布族。

具体关于指数型分布式的细节请看笔记概率统计随机过程之指数型分布族应用.md

自然指数族

我们在式（5）中给出了指数型分布族的一般形式 \[ p(x;\eta)=h(x)\exp\left\{\eta^T\mathbf{T}(x)-A(\eta)\right\}\tag{5} \] 但是对于广义线性模型的应用场景而言，还是复杂了一些，因此有一种简化的自然指数族。在自然指数族中，\(\mathbf{T}(\mathbf{x})=\mathbf{x}\)，不存在类似于\(x^2,x^3,\log(x),\frac{1}{x}\)这种带有函数关系的充分统计量，其可以简化写成： \[ p(x;\eta)=h(x)\exp\left\{\eta^T\mathbf{x}-A(\eta)\right\}\tag{6} \] 二项分布，负二项分布，伯努利分布，泊松分布，参数\(\alpha\)已知的Gamma分布，已知方差的高斯分布，参数\(\lambda\)已知的逆高斯分布（又称Wald分布）等都可以写成自然指数族形式，其他分布如卡方分布、Beta分布、帕累托分布，对数正态分布，一般正态分布，一般Gamma分布则无法写成自然指数族的形式。他们是否是自然指数族的核心就在于是不是充分统计量\(T(x)=x\)。

指数分散族

在自然指数族的基础上，研究者们为了方便探究分布的期望和方差，对自然指数族做了少些变形得到指数分散族。其处理方法是将自然指数族的规范形式(式(6))的规范（自然）参数\(\eta\)拆分成与位置（期望）相关的位置函数\(b(\vartheta)\)以及和方差相关的分散函数\(a(\phi)\)。其形式如下： \[ p(x;\vartheta)=\exp\{\frac{\vartheta^T x-b(\vartheta)}{a(\phi)}+c(x,\phi)\}\tag{7} \] 这种形式的指数族通常被称为指数分散族(exponential dispersion family,EDF)，\(a(ϕ)\)称为分散函数(dispersion function)，是已知的。\(ϕ\)称为分散参数(dispersion parameter)。\(\vartheta\)仍然叫自然参数(natural parameter)或者规范参数(canonical parameter)，它和自然指数族中参数差了个系数，因为两种模式中\(\vartheta^T x,\eta^Tx\)的模式都是参数\(\times\)充分统计量，所以不难发现，实际上我们对自然参数做一个\(\frac{1}{a(\phi)}\)倍的缩放。需要指出的是，在广义线性模型中，\(a(\phi)\)一般是 已知的，且通常是个常数系数，如果样本之间的重要性没有区别，我们可以令\(a(\phi)=\phi\)，即 \[ p(x;\vartheta)=\exp\{\frac{\vartheta^T x-b(\vartheta)}{\phi}+c(x,\phi)\}\tag{7.1} \]

指数分散族形式本质上是对自然指数族的参数\(\eta\)进行了拆分，把期望参数和方差参数拆分开（二者实际是可逆的变换）。使得自然参数\(\vartheta\)仅和期望\(μ\)相关，分散参数\(ϕ\)和分布的方差参数相关。分拆后，规范参数\(\vartheta\)仅和分布的期望参数\(μ\)相关，并且和\(μ\)之间存在一一映射的函数关系，换句话说，\(\vartheta\)和\(μ\)可以互相转化。 \[ \vartheta=f(\mu)\\ \mu=f^{−1}(\vartheta)\tag{8} \] 这一点接下来会证明。

由笔记概率统计随机过程之指数型分布族应用可知，指数分散族的期望和方差可表达为： \[ E[X]=b'(\vartheta)=\mu\tag{9} \] \[ \mathrm{Var}[X]=a(\phi)b''(\vartheta)\tag{10} \] 从期望和方差的关系，我们能发现\(\vartheta\)与\(\mu\)也是一一对应关系。根据式（9）可知，\(\vartheta\)与\(\mu\)有函数关系，且由于\(b'(\vartheta)\)的导数\(b''(\vartheta)\)是方差（恒大于0）乘以一个已知数\(a(\phi)\)（式（10）结论），因此\(b'(\vartheta)\)的导数必然恒为正数或负数（取决于已知数\(a(\phi)\)），即\(b'(\vartheta)\)必为单调函数，而单调函数必存在反函数，推得必存在\(b'^{-1}\)，使得\(\vartheta=b'^{-1}(\mu)\)。因此\(\vartheta\)与\(\mu\)是一一对应的。

我们定义配分函数\(b(\vartheta)\)的二阶导数为方差函数(variance function)，方差函数是一个关于期望\(μ\)的函数，即 \[ b''(\vartheta)=\nu(μ)\tag{11} \] 方差函数\(ν(μ)\)存在两种情况：

1. 方差函数是一个常量值，\(ν(μ)=b''(\vartheta)=C\)，此时分布的方差与均值无关。典型的分布就是正态分布。
2. 方差函数是一个关于均值\(μ\)的函数，\(ν(μ)=b''(\vartheta)\)，此时分布的方差与均值有关。

我们从一般的概率分布推导出指数族，然后其中的一个子集自然指数族，最后给它做一个变型的指数分散族，这么兜兜绕绕就是为了方便广义线性的计算与推导。

广义线性模型与指数型分布族（核心）

我们在前面提到了用概率的形式理解回归模型和指数型分布族，这两个概率论的概念其实是广义线性模型的核心，然而在实际应用中，不论是一般线性模型还是广义线性模型，都没有体现出概率的影子。概率论似乎从回归模型中消失了，这一节我们将概率论从幕后拖出来，展示其操作广义线性模型的真面目。

我们在从概率的角度解释回归章节中指出，根据参数估计形成的预测模型给出的预测值\(y_{pred}\)实际上并不会和观测到的\(y_{obs}\)(observation)完全一致。这其中的缘由主要有两点：

1. 观测值\(y_{obs}\)与对应的\(x_{obs}\)并不是严格的函数关系，而是在某个与\(x_{obs}\)相关的函数附近随机波动，即\(y_{obs}=f(x_{obs})+\varepsilon\)，其中\(\varepsilon\)是随机数，服从特定分布。但是，实际场景下\(f(x_{obs})\)与\(\varepsilon\)也不一定是相加的关系。
2. 我们的预测模型\(f(\cdot)\)并不保证一定准确，模型的参数也是通过观测数据通过统计推断的形式如（最大似然估计、最大后验概率估计）得到的，而非精确推导。

因此，我们认为响应变量\(y_{obs}\)也服从带有特定参数的分布\(Y=P(y_{obs}|x_{obs})\)，即\(y_{obs}\)是个随机变量。然而，\(y_{obs}\)不是一个独立的随机变量，它和\(x_{obs}\)有着密切关系，我们可以把\(x_{obs}\)看成\(Y\)的分布参数。实际应用中，我们不会说给出\(Y\)的分布当成\(y_{pred}\)让用户或者系统使用，这样一是不好用，二是计算起来经常无法得到数值结果。所以我们有个自然而然的想法：使用一个代表性的数字来替代\(Y\)的概率分布。

那么，如果只选一个数字来代替整个概率分布，大家的第一反应基本都是期望。因此，我们在建立回归模型的时候，希望根据给定的\(x_{obs}\)得到的预测值\(y_{pred}\)等于\(Y\)分布的条件期望\(E[Y|x_{obs}]\)。 \[ \mu=E[Y|x_{obs}]=y_{pred}=f(x_{obs})\tag{12} \] 式（12）是推导广义线性模型的核心。由于求概率的期望是一个抹除随机性的过程，因此在回归模型中，概率分布的痕迹被隐藏了起来。其中的\(f(\cdot)\)就是要训练的广义线性模型。

从函数的角度来看，\(x_{obs}\)为自变量，\(f(\cdot)\)为映射关系（广义线性模型），\(y_{pred}\)为因变量，又称响应变量（Response variable）。从式（12）我们也可以看出，建立的广义线性模型其根据自变量计算得到的因变量，应等于随机变量\(Y\)的期望。

我们要解析广义线性模型，首先就要从\(Y\)的分布谈起。经过之前章节的铺垫，我们应该已经猜到，\(Y\)是属于指数型分布族，即\(Y\sim P(y;\theta)\)。这是有现实意义的，比如认为估计的连续型随机变量属于高斯分布、二分类型随机变量属于伯努利分布等等。同时，我们也认识到，之所以广义线性模型都属于指数型分布族只不过是因为我们人为地挑了容易研究的这类分布族罢了。所以，从因果关系上来讲，并非广义线性模型都使用指数型分布族，而是我们先选中指数型分布族，然后把符合这些分布族的模型命名为广义线性模型。

虽说，\(P(y;\theta)\)是指数型分布族，但是广义线性模型用到的指数族模型并不像式（5）那么复杂，而是属于自然指数族或指数分散族。由于自然指数族和指数分布族是等效的，且指数分散族更适合模型的推导，同时我们一般不对样本的重要性有区分，因此我们下文主要使用\(a(\phi)\)为常函数的指数分散族的形式，如式（7.1）所示： \[ p(x;\vartheta)=\exp\{\frac{\vartheta^T x-b(\vartheta)}{\phi}+c(x,\phi)\}\tag{7.1} \] 为了区别观测数据\(x_{obs}\)和式（7.1）pdf中的自变量\(x\)，我们使用预测值\(y\)替代式（7.1）中的\(x\)，即 \[ p(y;\vartheta)=\exp\{\frac{\vartheta^T y-b(\vartheta)}{\phi}+c(y,\phi)\}\tag{7.2} \] 此式（7.2）即为在下文中使用的广义线性模型概率分布表达式。

根据公式（9）我们可知，给定\(\vartheta\)的\(Y\)的期望为配分函数\(b(\vartheta)\)的一阶导数: \[ E[Y;\vartheta]=b'(\vartheta)=\mu\tag{9.1} \] 也就是说，我们使用期望\(b'(\vartheta)\)来代表\(Y\)的整个分布，并且也应是广义线性模型预测的结果\(y_{pred}\)。需要指出的是，\(\vartheta\)是\(x_{obs}\)的函数（有确定关系），因此\(E[Y|x_{obs}]=E[Y|\vartheta]\)。综合式（12）（9.1）我们有： \[ b'(\vartheta)=\mu=y_{pred}=f(x_{obs})\\ \Rightarrow b'(\vartheta)=f(x_{obs})\tag{13} \] 从式（12）、式（9.1）、式（13）我们通过求期望的方式去除了\(Y\)的随机性，得到了一个确定性的表达式，这也将概率统计的影子从广义线性模型中消去了。如果说还是留有概率分布的痕迹的话，那么就只有指数分散族的配分函数的导数\(b'(\vartheta)\)还在其中。

式（13）还说明我们要求的广义线性模型的表达式和配分函数的导数\(b'(\vartheta)\)是存在密切关系的！ 此外，我们在式（9）（10）的介绍中指出由于\(\vartheta\)与\(\mu\)是一一对应的，即存在反函数 \[\vartheta=b'^{-1}(\mu)\tag{14}\] 这个式子将会把广义线性模型中的“线性”和指数型分布族联系起来。

线性的体现

我们前面说了指数型分布族、期望，甚至得到了自变量\(x_{obs}\)与\(y_{pred}\)的某种关系\(\mu=b'(\vartheta)=y_{pred}=f(x_{obs})\)，然而还有一个关键点没有解决，那就是如何将概率分布的参数\(\vartheta\)和\(x\)观测值\(x_{obs}\)联系起来。我们可以通过上式发现，\(\vartheta\)和\(x_{obs}\)是通过均值\(\mu\)联系在一起的。而\(\vartheta\)与\(\mu\)的关系可以通过式（14）确定。然而，\(x_{obs}\)与\(\mu\)的关系还没有定下来。这也是广义“线性”模型中线性一词的由来，我们人为地设计一个线性预测器： \[ \kappa=\beta^T x_{obs}+b\tag{15} \] 即我们设定参数\(\kappa\)是\(x_{obs}\)的线性组合，\(x_{obs}\)可以是标量也可以是向量。为了简洁性，通常会人为的为\(x\)扩充一个一维常量值1，并且把截距参数\(b\)算在\(β\)中，这样上述线性函数可以写成向量內积的形式。 \[ \kappa=\beta^T x_{obs}\tag{15.1} \] 结合式（13）（14）我们可以得到 \[ \mu=b'(\vartheta)=h(\kappa)\tag{16} \] \[ \vartheta=b'^{-1}(\mu)=b'^{-1}(h(\kappa))=b'^{-1}(h(\beta^T x_{obs}))\tag{17} \] 根据式（13）的关系，显然有\(f(x_{obs})=h(\kappa)=h(\beta^T x_{obs})\)。由于\(f\)与\(h\)之间自变量\(x_{obs},\beta^T x_{obs}\)是线性变换，因此只要\(f(\cdot)\)存在，那么一般情况下，\(h(\cdot)\)必然存在。这样我们就得到了广义线性模型，我们使用式（16）来通过观测值\(x_{obs}\)与\(\beta\)的线性组合得到预测值\(y_{pred}=\mu\)。而式（17）则表明了线性预测器\(\beta^T x_{obs}\)与自然参数\(\vartheta\)的关系。现在还有一个问题，就是函数关系\(h(\cdot)\)是什么样的呢？

连接函数与激活函数

从式（16）（17）我们可以看出预测值（期望）和线性预测器之间的关系。我们将令 \(g(\mu)=h^{-1}(\mu)\)称为连接函数（Link function），之所以叫这个名字是因为它将线性预测器和最终的预测值（期望）联系到了一起，即 \[ g(\mu)=h^{-1}\circ h(\beta^T x_{obs})=\beta^T x_{obs}=\kappa\tag{18} \] 但是在实际使用中，连接函数并不常用，反倒是连接函数的反函数\(g^{-1}=h\)更常见，因为我们可以用过\(g^{-1}(\beta^T x_{obs})\)计算出广义线性模型的预测值\(y_{pred}\)，我们称\(h=g^{-1}\)为激活函数(Activation function)。 \[ y_{pred}=g^{-1}(\kappa)=g^{-1}(\beta^T x_{obs})=h(\beta^ T x_{obs})\tag{19} \] 显然，广义线性模型中通过线性预测器和激活函数就可以得到预测结果\(y_{pred}\)，而激活函数与连接函数是反函数的关系。现在，如果我们能知道连接函数，就可以完成整个广义线性模型了！

那么现在就有一个问题，我们如何来确定连接函数呢？此外，也引出另一个问题:连接函数是唯一的吗？我们先回答后一个问题，即连接函数并不是唯一的。下面我们讨论如何找连接函数\(g\)。

前面我们提到过\(\vartheta\)是\(x_{obs}\)的函数（有确定关系），且有\(E[Y|x_{obs}]=E[Y|\vartheta]=\mu=h(x_{obs})\)，因此本质上期望\(\mu\)也是一个关于\(x_{obs}\)的函数。那么，式（18）可简写成 \[ g(\mu)=\beta^T x_{obs}\tag{18.1} \] 这是一个比较有意思的式子，左右两边根源自变量都是\(x_{obs}\)。那么连接函数\(g\)就必须要满足2个条件：

1. 由于\(\beta^T x_{obs}\)的取值可能是整个实数域，但是期望\(\mu\)是有范围的，比如泊松分布的均值必大于0，即\(\mu\in (0,+\infty)\)，因此\(g\)必须将\(\mu\)的取值范围映射到整个实数域，这称为定义域要求。连接函数本质上，就是把实数域范围的\(\beta^T x_{obs}\)转换到特定分布合法的\(\mu\)值空间上。
2. 此外，我们希望\(g\)是可逆的（可微且严格单调），这样\(\mu\)与\(\beta^T x_{obs}\)就有了一一对应的关系，这是因为函数一旦不可逆\(\beta^T x_{obs}\)就可能求出多个预测值，这显然是不符合实际情况的，这称为可逆要求。

满足以上2点，就可以初步得到一个连接函数。可以看出，上面2个要求还是相对宽泛的，不同的映射可能得到不同的线性关系组合即\(\beta_1^Tx_{obs},\beta_2^Tx_{obs},\dotsb\)，这些都满足连接函数的要求，也就是说满足定义域映射关系的可逆函数都可以作为连接函数，因此连接函数并不是唯一的。

举两个例子：

一：伯努利分布，我们知道其期望\(\mu\in(0,1)\)，那么连接函数\(g\)应该能够将定义域为\((0,1)\)的范围映射到整个实数域\(\R\)，那么下面三个连续单调函数（可逆）都可以满足：

1. Logit：\(\log(\frac{\mu}{1-\mu})\)；
2. Probit：\(\varPhi^{-1}(\mu)\)，其中\(\varPhi^{-1}\)是正态分布概率累计函数的反函数，实际上在伯努利分布中，任意定义域为\(\R\)的概率累计函数的反函数都可以作为连接函数；
3. 互补Log-Log：\(\log(-\log(1-\mu))\)

可以验证以上三个函数都可以将\(\mu\in(0,1)\)映射到实数域\(\R\)。

二：泊松分布，泊松分布的期望值\(\mu>0\)，那么连接函数\(g\)需要将定义域\((0,\infty)\)的范围映射到整个实数域\(\R\)，那么下面两个连续单调函数（可逆）都可以满足：

1. Log函数：\(\log(\mu)\);
2. 初等函数：\(\mu-\frac{1}{\mu}\)，虽然此函数在全局并不是连续单调的，但是在\(\mu>0\)时是全局单调的，存在反函数\(y=\frac{x+\sqrt{x^2+4}}{2}\)。

既然连接函数并不是唯一的，那么我们就要从中选出最为合适的连接函数。下一节，我们介绍一个比较常用的连接函数选择方式：标准连接函数。

标准连接

我们之前说个，连接函数并不是唯一的，那么选择哪一个连接函数比较好呢？这里我们推荐一种比较自然的选择，注意只是说比较自然，但不一定是最合适的，在不同场景下某些连接函数确实表现地比其他连接函数更好。

所谓比较自然的选择，就是选择连接函数使得\(g(\mu)=\vartheta=\beta^Tx_{obs}\)。首先，\(\vartheta\)的取值范围就是\(\R\)，因此可以当成连接函数的值域，同时，我们也能保证常见情景下，这样选择的连接函数是可逆的（可用指数型分布族配分函数的二阶导数含义证明）。选择标准连接函数的一大优势就是可以极大地简化数学运算，非常契合广义线性模型的最大似然估计，我们在下一章节会详细讨论这个问题，这里先记住这个结论。

根据式（16）与标准连接函数的条件\(g(\mu)=\vartheta=\beta^Tx_{obs}\)可以推出： \[ g(\mu)=g(b'(\vartheta))=\vartheta\tag{20} \] 其中，\(\vartheta\)经过两次变换\(g\circ b'\)又变回了\(\vartheta\)，显然有\(g\circ b'\)为恒等变换，那么\(g\)与\(b'\)就互为反函数！即有 \[ g^{-1}=b'\tag{21} \] 这样我们的标准连接函数就可以根据指数型分布族的配分函数\(b'(\vartheta)\)的反函数推得了！我们终于不用再在一堆连接函数中盲目地搜索、构造，只需要通过\(b'(\vartheta)\)就能求得。同时，我们知道，激活函数与连接函数也是互为反函数，那么根据式（21）的关系也不难发现： \[ b'=h\tag{22} \] 指数型分布族的配分函数的一阶导数正好就是激活函数的形式！我们最终得到了一个既简洁又优美的结果！

下表列出了常用广义线性模型的配分函数、激活函数、连接函数以及典型使用场景。注：在自然指数族中，我们只允许充分统计量\(T(x)=x\)，在其他文献中有别的设置方式，因此标准连接函数会略有所不同，大多数只是系数的差别。

广义线性模型的参数关系

我们在上表中，可以看出广义线性模型的表达式等于激活函数，即\(f_{GLM}(x)=g^{-1}(\beta^T x)\)，如果我们梳理一下上面的内容，可总结出广义线性模型的三大组成部分。

1. 一个指数族分布(指数分散族)作为响应变量\(Y\)概率分布\(p(Y;\vartheta)\)，被称为随机组件(random component)。
2. 一个线性预测器\(\vartheta=β^Tx_{obs}\)，被称为系统组件(systematic component)。
3. 一个连接函数\(g\)使得\(E[Y]=\mu=g(\vartheta)=g(\beta^T x_{obs})\)，描述了随机组件和系统组件之间的关系。

用图像表达他们的关系为（假设使用标准连接函数）：

广义线性分布关系图

先从图坐往右看，待定参数\(\beta\)与观测值\(x_{obs}\)组成线性预测期，其结果为指数分散族的参数\(\vartheta\)，且其分散参数\(\phi\)已知。由于指数分散族和常使用的一般概率分布表达式是可以相互转换的，我们也可以得到原参数\(\theta\)和自然参数\(\vartheta\)的转换关系，可以证明它们是一一对应的。

再从图右往左看，由于观测结果\(y_{obs}\)与\(x_{obs}\)存在相关性，而非确定函数关系，我们认为其\(y_{obs}\)服从特定的分布，我们认为这种分布是指数型分布族（指数分散族），因此可以用\(P(y;\vartheta)\)来表达\(y_{obs}\)的分布，用随机变量\(Y\)表示。当然，这个随机变量也可以用原始参数的概率分布\(P(y;\theta)\)表示。但是，在实际使用中，我们更希望得到一个具体值，而非一个分布，因此我们应使用最具代表性的值——期望\(\mu\)来表示分布，即期望应等于广义线性模型的预测值\(\mu=y_{pred}\)。同时，指数分散族有一个非常好的性质，其配分函数的导数\(b'(\vartheta)\)正好等于期望\(\mu\)。而原始参数的概率分布\(P(y;\theta)\)与期望\(\mu\)的关系\(\varphi(\theta)\)则需要根据具体情况计算，因此在广义线性模型中不太常用。

由于\(\mu\)和线性预测期直接不是相等关系，需要一个桥梁，即连接函数\(g(\mu)\)，它阐述了期望和线性预测期之间的关系：\(g(\mu)=\vartheta=\beta^T x_{obs}\)。我们可以证明这个连接函数必存在反函数\(g^{-1}\)，这样我们就可以使用这个反函数来预测结果，即\(g^{-1}(\vartheta)=g^{-1}(\beta^T x)\)。我们将这个反函数称之为激活函数。最后我们发现配分函数的导数\(b'(\vartheta)\)与激活函数\(g^{-1}(\vartheta)\)有着一样的作用，因此我们得到激活函数就是配分函数的导数，即\(b'=g^{-1}\)。它们就是广义线性模型的函数\(f(x_{obs})=b'(\beta^T x_{obs})=g^{-1}(\beta^T x_{obs})\)。

广义线性模型的最大似然估计

本笔记的最后分两步，首先我们将说明使用标准连接函数时，使用最大似然估计方法我们可以得到非常简洁的解析结果；然后我们将首位呼应，回收文章开头引子中提出的疑问。

标准连接函数下的广义线性模型最大似然估计

最大似然估计是广义线性模型中参数估计的一般方法，其目的是最大化似然函数，具体细节可参考笔记《概率统计随机过程之最大似然估计拓展.md》。其模型通常假设如下：我们有一组\(n\)个观测值\((x_i,y_i)\in(\R^p\times \R),i=1,2,\dotsb,n\)，其中\(x_i\)是一个\(p\)维向量，\(y_i\)是观测结果是一个实数。观测值\(x_i\)会影响到广义线性模型的参数\(\vartheta\)，或者说\(\vartheta\)是\(x_i\)的函数，即\(\vartheta=\vartheta(x_i)\)。在最大似然估计中，我们假设各个观测值之间是独立的，那么每一个观测值发生的概率即为： \[ p(y_i)=\exp\{\frac{\vartheta(x_i)y_i-b(\vartheta(x_i))}{\phi}+c(y_i,\phi)\}\tag{23} \] 显然，观测值\(x_i\)通过其因变量\(\vartheta(x_i)\)影响概率与期望，\(x_i\)给定时\(\vartheta\)也确定，然后\(b(\vartheta)\)及其一阶导数也确定了，最终可以确定\(\mu=b'(\vartheta)\)。我们先将独立的观测值相乘取\(log\)得到对数似然函数： \[ \begin{aligned} l_n(\vartheta(X);Y)&=\log \prod_{i=1}^n \exp\{\frac{\vartheta(x_i)y_i-b(\vartheta(x_i))}{\phi}+c(y_i,\phi)\}\\ &=\sum_{i=1}^n \frac{\vartheta(x_i)y_i-b(\vartheta(x_i))}{\phi}+c(y_i,\phi) \end{aligned}\tag{24} \] 下面我们根据参数\(\vartheta\)，期望\(\mu\)、线性预测器\(\beta^T x_i\)以及连接函数\(g\)之间的关系对式（24）进行变形： \[ \left .\begin{aligned} 式(9.1)\Rightarrow \vartheta=b'^{-1}(\mu)\\ 式(18.1)\Rightarrow \mu=g^{-1}(\beta^T x_i) \end{aligned}\right\}\Rightarrow \vartheta(x_i) = b'^{-1}\circ g^{-1}(\beta^T x_i)\\ \Rightarrow \vartheta(x_i)=(g\circ b')^{-1}(\beta^T x_i)\tag{25} \] 上式集中体现了\(\vartheta\)与\(x_i\)之间的关系，代入式（24）可将变量\(\vartheta(x_i)\)替换成\(\beta\)： \[ l_n(\beta;Y;X)=\sum_{i=1}^n \frac{(g\circ b')^{-1}(\beta^T x_i)y_i-b[(g\circ b')^{-1}(\beta^T x_i)]}{\phi}+c(y_i,\phi)\tag{26} \] 在场景中，\((x_i,y_i)\)都是已经获得的观测值，只有参数\(\beta\)是未知的，最大似然估计就是求\(\beta\)使得对数似然函数（等效于似然函数）最大。

一般情况下，式（26）的计算并不是平凡的，但是如果我们选用标准连接函数，这个式子就能得到极大简化！根据式（21）有：\(g=b'^{-1}\)，这使得式（26）中的函数复合\(g\circ b'=I\)，即二者抵消，是一个恒等变换！此时，式（26）可简化为： \[ l_n(\beta;Y;X)=\sum_{i=1}^n \frac{\beta^T x_iy_i-b(\beta^T x_i)}{\phi}+c(y_i,\phi)\tag{26.1} \] 其中，\(\beta^T x_iy_i\)是一个简单的线性函数；由式（10）可知\(b(\beta^T x_i)\)二阶导海森矩阵正定是一个严格凸函数；而最后的\(c(y_i,\phi)\)与\(\beta\)无关，求导时不造成影响。这使得采取标准连接函数的最大似然估计可以通过凸优化的方法求得唯一极值，大大简化优化流程，免除了函数复合\(g\circ b'\)在求导计算中复杂情形，这也是我们愿意选用标准连接函数的核心原因。

附标准连接函数下最大似然估计的一二阶导数：

一阶导数： \[ \nabla_\beta l_n(\beta;Y;X)=\sum_{i=1}^n \frac{y_ix_i-b'(\beta^T x_i)x_i}{\phi};x_i\in \R^p\tag{27} \] 二阶导数： \[ \nabla^2_\beta l_n(\beta;Y;X)=\sum_{i=1}^n\frac{-b''(\beta^T x_i)x_ix_i^T}{\phi};x_i\in\R^p;x_i x_i^T \text{positive definite}\tag{28}\\ \forall y\neq 0\in \R^p\quad y^T x_i x_i^T y=(x_i^T y)^T (x_i^Ty)>0 \]

Logistics回归优化举例

伯努利分布又叫两点分布或者0-1分布，是最简单的概率分布形式之一，其对应的广义线性模型就是Logistics回归，即二项分类。常见伯努利分布写成： \[ p(y;\theta)=\theta^y(1-\theta)^{1-y},y\in\{0,1\},\theta\in[0,1] \] 转写为指数型分布族形式为： \[ \begin{aligned} p(y;\theta)&=\exp\{y\ln{\theta}+(1-y)\ln{(1-\theta)}\}\\ &=\exp\{y\ln(\frac{\theta}{1-\theta})+\ln{(1-\theta)}\} \end{aligned} \] 令\(\vartheta=\ln(\frac{\theta}{1-\theta})\)将其转换成指数分散族： \[ p(y;\vartheta)=\exp\{y\vartheta-\log(1+e^\vartheta)\} \] 其\(n\)个观测值组成的对数似然函数为： \[ l_n(\vartheta(X);Y)=\sum_{i=1}^n (y_i\vartheta(x_i)-\log(1+e^{\vartheta(x_i)})) \] 我们采用标准连接函数： \[ \vartheta(x_i)=\beta^T x_i \] 则对数似然函数可写为： \[ l_n(\beta;Y;X)=\sum_{i=1}^n (\beta^T x_iy_i-\log(1+e^{\beta^Tx_i})) \] 对\(\beta\)求一阶，二阶导数分别为： \[ \nabla_\beta l_n(\beta;Y;X)=\sum_{i=1}^n (y_ix_i-\frac{e^{\beta^Tx_i}}{1+e^{\beta^T x_i}}x_i);x_i\in \R^p\\ H_{l_n}(\beta)=\nabla^2_\beta l_n(\beta;Y;X)=\sum_{i=1}^n \frac{e^{\beta^Tx_i}}{(1+e^{\beta^T x_i})^2}x_ix_i^T;x_i\in\R^p; \] 注意，\(x_ix_i^T\)是由向量张成的矩阵。那么根据牛顿法，其优化迭代步骤为： \[ \beta^{(k+1)}=\beta^{(k)}-H_{l_n}(\beta^{(k)})^{-1}\nabla_\beta l_n(\beta^{(k)};Y;X) \]

补充：牛顿法的简化方法之一Fisher分数法

由于计算海森矩阵求和这一步很繁琐，因此我们可以用期望来替代求和操作（具体原理参见笔记概率统计随机过程之最大似然估计拓展中最大似然估计与相对熵（KL散度）、交叉熵的等价性那一节）。因此有： \[ H_{l_n}(\beta)=\nabla^2_\beta l_n(\beta;Y;X)=\sum_{i=1}^n\frac{-b''(\beta^T x_i)x_ix_i^T}{\phi}\\ \simeq E[H_{l_n}(\beta)]=-I(\beta) \] 其中，\(I(\beta)\)是Fisher信息量。所以海森矩阵也可以被Fisher信息矩阵替代，即为： \[ \beta^{(k+1)}=\beta^{(k)}+I(\beta^{(k)})^{-1}\nabla_\beta l_n(\beta^{(k)};Y;X) \] 这种方法称为Fisher分数法，这算是一种近似的拟牛顿法，其收敛速度和牛顿法相近，但是计算量降低。

回答引子的疑问，最大似然估计形势下的迭代优化

还记得在文章开头时的引子吗？

如果刚学完线性回归和Logistics回归，那么是否会注意到，二者的梯度更新步骤都是(虽然\(h_{\vec\theta}(\vec x^{(i)})\)的定义不同)： \[ \theta_j=\theta_j-\alpha(h_{\vec\theta}(\vec x^{(i)})-y^{(i)})x_j^{(i)}\\ h_{\vec\theta}(\vec x^{(i)})=\begin{cases} \vec{\theta}^T \vec{x},\quad线性回归\\ \frac{1}{1+e^{-\vec{\theta}^T \vec{x}}},\quad Logistics回归\end{cases} \] 其中，\(\vec\theta, \vec x^{(i)}\)分别是参数向量，第\(i\)个观测数据的向量。下标\(j\)表示第\(j\)个分量，\(\alpha\)表示更新的步长。

这一是因为我们通过指数型分布族（指数分散族）给出了广义线性模型响应变量\(Y\)的统一形式，第二点就是由于我们使用的是最大似然估计，它的似然函数求导步骤与指数分散族非常契合。实际上引子中的例子就是使用标准连接函数时一阶导数的应用。在梯度下降法中，通用的迭代公式为： \[ \theta_j=\theta_j-\alpha\frac{\partial }{\partial \theta_j}L(\theta) \] 代入式（27）即为(变量\(\beta\)替换为\(\theta\)，另外最大似然估计要变成最小化损失函数，加个负号): \[ \theta=\theta-\alpha\sum_{i=1}^n -\frac{y_ix_i-b'(\theta^T x_i)x_i}{\phi}=\theta-\alpha\sum_{i=1}^n(b'(\theta^T x_i)-y_i)x_i/\phi \] 其中，\(b'(\theta^T x_i)\)写成激活函数形式为\(h_{\vec\theta}(\vec x^{(i)})\)，采用单步迭代不需要求和\(\sum\)，如果分散系数为1，那么最终结果的各维度分量就是引子中的形式。

广义线性模型的求解（IRLS算法）

通常情况下，如果我们不考虑计算的复杂程度，广义线性模型可以使用梯度下降法或者牛顿进行求解。在这里，我们介绍一种广义线性模型的较常用的优化算法Iteratively Re-wighted Least Squares算法。这个算法能够将循环求和运算变换成矩阵运算，提供一个较为整体化的便捷计算形式。

TODO

参考文献

1. https://zhangzhenhu.github.io/blog/glm/source/%E5%B9%BF%E4%B9%89%E7%BA%BF%E6%80%A7%E6%A8%A1%E5%9E%8B/content.html#id2
2. https://www.youtube.com/watch?v=mc1y8m9-hOM&list=PLUl4u3cNGP60uVBMaoNERc6knT_MgPKS0&index=21