8. 参数估计¶
本章我们讨论 GLM 模型的参数估计算法，
我们将统一的以指数族的形式展现算法过程，这样适用于指数族中的所有具体分布，
我们的目标是让大家对 GLM 的基础理论有个全面的了解，同时我们会着重强调算法成立的假设及其一些限制。
传统上，对于单参数的指数族分布可以运用梯度下降法和牛顿法进行参数估计，
梯度下降法的优点是算法实现简单，缺点是收敛速度不如牛顿法。
梯度下降法和牛顿法在形式上是非常相似的，二者都是沿着目标函数的负梯度方向寻找最优解，
不同的是传统梯度下降法利用一阶导数，而牛顿法利用二阶导数，牛顿法相对于梯度下降法收敛速度会更快，
但是由于二阶导数的引入也使得牛顿法的计算复杂度增加很对，甚至很多时候无法计算。
在用牛顿法对指数族模型进行参数估计时，不同的分布拥有不同的梯度表达式，所以每种分布都需实现一个适合自己的牛顿法。
这里我们同时会介绍牛顿法的一个变种算法，迭代重加权最小平方法(iteratively reweighted least squares,IRLS)。
GLM 框架下的模型都可以以统计的形式运用 IRLS 算法进行参数估计，这是 GLM 非常有吸引力的一点。
IRLS 算法的另一个特点是不需要对估计参数 \(\hat{\beta}\) 进行初始化设置。
8.1. 最大似然估计¶
最大似然估计是应用十分广泛的一种参数估计方法，
其核心思想是通过极大化最大似然函数找到参数的最优解，
GLM 中参数估计就是使用的最大似然方法。
注意在 GLM 中，最大似然方法不能同时估计线性协变量参数 \(\beta\) 和分散参数
\(\phi\) ，
在GLM中的最大似然估计通常是假设分散参数 \(\phi\) 已知的情况下，
估计协变量参数 \(\beta\) 。
假设响应变量 \(Y\) 是 GLM 中的指数族分布，
协变量 \(X\) 和响应变量 \(Y\) 的一个样本集为 \(\mathcal{D}=\{(x_1,y_1),(x_2,y_2),\dots,(x_N,y_N)\}\)
，样本之间是相互独立的，本书中的最大似然估计都是建立在样本独立的假设之上。
所有样本的联合概率为
\(\theta\) 是指数族分布的自然参数，\(\phi\) 是分散参数。 样本的联合概率又被称为样本集的似然函数，
公式(8.1.1) 和 公式(8.1.2) 的区别在于， 前者是一个概率密度(质量)函数，是在给定 \(\theta,\phi\) 的条件下关变量 \(Y\) 的函数； 后者是一个似然函数，其表达是在给定观测样本 \(y\) 的条件下关于未知参数 \(\theta,\phi\) 的函数。
因为似然函数是一个连乘形式，所以通常我们会对其进行一个对数转换(log-transform)进而得到一个连加的形式， 连加的形式更方便进行计算。 连乘变成连加有两个好处，(1)更容易求导和极大化操作；(2)似然函数是概率连乘，而概率都是小于1的，大量小于1的数字连乘产生更小的数字， 甚至趋近于0，而计算机的浮点数精度是通常无法处理这么小的数字的，所以加对数更方便计算机进行数值处理。
加了对数的似然函数被称为对数似然函数，通常用符号 \(\ell\)
表示，对数似然函数是最大似然估计(ML)的核心，
GLM 模型的对数似然估计函数为
\(\theta_i\) 是自然参数， \(b(\theta_i)\) 是累积函数，它描述的是分布的矩(moment)； \(\phi\) 是分散参数(dispersion parameter)，影响着分布的方差； \(c(\cdot)\) 是归一化项。 归一化项不是 \(\theta\) 的函数，而是简单地缩放基础密度函数的范围使得整个函数的积分（或求和）为1。 在上一节我们讨论过，指数族分布的期望与方差可以通过 \(b(\theta)\) 的导数求得。
并且我们知道，均值参数 \(\mu\) 和自然参数 \(\theta\) 是存在一个可逆的函数关系的，也就是说 \(\mu\) 可以看做是关于 \(\theta\) 的一个函数，反之，\(\theta\) 也可看做是一个关于 \(\mu\) 的函数。 基于这个事实，我们可以把 \(b''(\theta)\) 看做是一个关于 \(\mu\) 的函数，记作 \(\nu(\mu)\) 。
因此，方差 \(V(Y)\) 就可以被看成是函数 \(\nu(\mu)\) 和分散函数 \(a(\phi)\) 的乘积，通常我们把 \(\nu(\mu)=b''(\theta)\) 称为方差函数(variance function)， 注意：虽然叫方差函数，但方差函数的值不是方差本身。 有时 \(b''(\theta)\) 会是一个常数量(constant)，比如高斯分布，此时分布的方差为：
这时，分布的方差就不会受到均值的影响了。 另外方差函数 \(\nu(\mu)\) 可以通过简单方式求得。
显然，当 \(\mu\) 与 \(\theta\) 之间的映射函数是线性函数时，一阶偏导 \(\frac{\partial \mu}{\partial \theta}\) 就是一个常数值。另外，我们知道反函数的导数就等于原函数导数的倒数，所以有：
在 GLM 框架下，输入变量 \(X\) 和其系数 \(\beta\)
组成一个线性预测器 \(\eta=\beta^Tx\) 。
\(\eta\) 和分布的均值(期望)通过连接函数(已知的)连接在一起 。
其中 \(\beta\) 和 \(x\) 都是一个向量，\(\beta^Tx = \sum_j \beta_j x_j\) ， 因此有：
线性预测器 \(\eta\) 的值空间并没有特别的限定， 其值空间是整个实数域 \(\eta \in R\) 。 而 \(\mu\) 的取值范围是特定分布相关的， 不同指数族分布，\(\mu\) 的取值范围是不同的，比如高斯分布 \(\mu \in R\) ，二项分布 \(\mu \in [0,1]\) 。 因此，连接函数的一个目的就是将线性预测器的值映射到响应变量期望参数的范围。
现在让我们回到最大似然估计，最大似然估计的思想是使得似然函数取得最大值的参数值为模型的最优解。 根据微分理论，函数取得极值的点其一阶偏导数为 \(0\)， 然而导数为 \(0\) 的点不一定是最大值的点，也可能是驻点、最小值点， 所以最大似然估计要求似然函数最好是凹函数。 利用求导的链式法则对GLM模型的对数似然函数进行求导， 注意，参数 \(\beta\) 是一个向量，所以这里是偏导数，
其中 \(i\) 是观测样本的编号，\(j\) 是参数向量的下标。 \(x_{ij}\) 表示第 \(i\) 条观测样本的第 \(j\) 列特征值， \(y_i\) 是响应变量的观测值，\(a(\phi)\) 通常认为是已知的。 \(\mu_i\) 是 \(y_i\) 的期望，也是模型的预测值， 方差函数 \(\nu(\mu_i)\) 是关于 \(\mu_i\) 的函数，因此也可以算出。 \(\frac{\partial \mu}{\partial \eta}\) 是响应函数(或者说是激活函数)关于 \(\eta_i\) 的导数， 在确定了连接函数的形式后也是可以算出的。
公式(8.1.11) 是 GLM 标准形式下对数似然函数的一阶偏导数，
GLM 框架下的任意模型都可以按照此公式计算偏导数，
只需要按照特定的分布和连接函数替换相应组件即可。
对数似然函数的一阶导数又叫得分统计量(score statistic,Fisher score)，或者得分函数(score function)， 常用符号 \(U\) 表示。
\(U\) 的表达式中只有 \(y_i\) 是随机变量的样本，其它都是数值变量， \(U\) 是一个关于样本的函数，所以它是一个统计量(statistic)， 得分函数(score function)有时也叫作得分统计量(score statistic)， 统计量也是随机变量。 我们知道 \(\mathbb{E}[y_i]=\mu_i\)， 而统计量 \(U\) 是变量 \(y\) 的函数，因此 \(U\) 期望值为：
统计量 \(U\) 的方差 \(\mathcal{J}=\mathbb{E}\{(U-\mathbb{E}[U])(U-\mathbb{E}[U])^T \}=\mathbb{E}[UU^T]\) 又被称为费希尔信息(Fisher information)，或者信息矩阵(information matrix)。
\(\mathbb{E}[(y_i-\mu_i)(y_l-\mu_l)]\) 是 \(y_i\) 与 \(y_l\) 的协方差， 根据样本独立性假设，有 \(y_i \perp \!\!\! \perp y_l (\ l \ne i)\) 成立，因此 \(y_i\) 与 \(y_l\) 的协方差为0， 即 \(\mathbb{E}[(y_i-\mu_i)(y_l-\mu_l)]=0\) 。而 \(\mathbb{E}_{y}[(y_i-\mu_i)^2]\) 表示 \(y_i\) 的方差，有 \(\mathbb{E}_{y}[(y_i-\mu_i)^2]=V(y_i)=a(\phi)\nu(\mu_i)\) 。最终化简为
在最大似然估计的理论中， 通过令 \(U=0\) 求得参数估计值，这个等式被称为估计等式(estimating equation)， 有的资料中也叫正规方程(normal equation)。
在这个方程中，有 \(\mu_i=r(\eta_i)=r(x_i^T \beta)\) ，函数 \(r(\dot)\) 是连接函数的反函数，称为响应函数，是已知的。 协变量系数 \(\beta\) 是方程的未知量，也是模型的未知参数， 是我们想要求解的。 分散函数 \(a(\phi)\) 通常被认为是已知的， 假设 \(a(\phi)=\phi\) ， 并且 \(\phi\) 与样本无关，即所有样本具有相同的值。 当 \(U_j=0\) 时，有
显然，在样本具有相同分散参数 \(\phi\) 的假设之下， 协变量参数 \(\beta\) 的最大似然估计是不受 \(\phi\) 影响的。
现在我们以传统线下回归模型为例，演示下如何利用估计等式进行参数求解。
传统线性回归模型也是 GLM 的一员，
相当于响应变量 \(y_i\) 是高斯变量，
\(y_i \sim \mathcal{N}(\mu_i,\sigma^2=1)\)
，并且连接函数是恒等函数 \(\eta_i=\mu_i\)
，响应函数作为连接函数的的反函数，自然也是恒等函数，即 \(\mu_i=\eta_i\)
，因此响应函数的导数是常量 \(1\)
。
传统线性回归模型中，方差是常量，\(V(y_i)=\sigma^2=1\) ，因此有
各项代入到估计方程中， \(U_j\) 简化为
上式是单个参数 \(\beta_j\) 的得分统计量 \(U_j\) ，转成向量为
移项可得参数的估计值为
我们发现标准连接函数的高斯模型，估计等式 \(U=0\) 可以得到解析解， 这是高斯模型独有的特性，其它模型或者连接函数是不具备这个特性的。
在 GLM 中，估计等式要想得到解析解，需要满足两个条件：
连接函数的是标准连接函数的。
连接函数是线性函数。
根据标准连接函数的定义，标准连接函数的使得 \(\theta_i=\eta_i\) ，此时有 \(\partial \theta_i / \partial \mu_i = \partial \eta_i / \partial \mu_i\) ，得分统计量 \(U\) 可以得到简化。
当采用标准连接函数时，得分统计量中，响应函数的导数和连接函数的的导数互相抵消掉， 得分统计量 \(U\) 得到简化。 再根据上文所述，\(a(\phi)\) 不影响参数估计结果，可以去掉。
转换成矩阵的形式为
移项可得
当 \(\mu\) 与 \(\eta\) 是线性关系时，比如 \(\mu_i = \alpha \eta_i = \alpha (x_i^T \beta)\) ， 公式(8.1.26) 才能求得解析解。
在 GLM 中，能同时满足这两个条件的，只有高斯模型，其它的模型都不符合第二点。
对于无法取得解析解的模型，可以用数值法求解。最常用的数值法有梯度下降法和牛顿法，
梯度下降法仅利用似然函数的一阶导数，而牛顿法同时利用似然函数的一阶导数和二阶导数，
下节我们介绍最大似然估计的数值求解法。
8.2. 泰勒级数¶
最大似然的求解需要求解正规方程，
然而在 GLM 中，正规方程并不是一定存在解析解的，需要满足一些限制条件才行，
解析解的方式并不具备通用性，我们需要采用更一般的方法，逼近法，也叫迭代法、数值法。
迭代法又可以简单分为一阶导(梯度下降法系列)和二阶导(牛顿法系列)，实际这两种都可以通过泰勒级数(Taylor series)进行推导。
泰勒级数有很多个名字，泰勒公式(Taylor formula)、泰勒级数(Taylor series)、
泰勒展开(Taylor explanation)、泰勒定理(Taylor theory)等，都是一回事。
设 \(n\) 是一个正整数。如果定义在一个包含 \(x_0\) 的区间上的函数 \(f\) ，在 \(x_0\) 处 \(n+1\) 次可导，那么对于这个区间上的任意 \(x\) ，都有
其中 \(f^{(n)}\) 表示 \(f\) 的 \(n\) 阶导数， 泰勒展开表达的就是 \(f(x)\) 可以用其附近的点 \(f(x_0)\) 近似的表示。
提示
注意，本书讨论的迭代求解算法默认目标函数都是凸函数，也就是函数有唯一的极值点。 关于非凸函数以及带约束的优化问题，请读者参考其它资料。
8.3. 梯度下降法¶
我们把对数似然函数按照泰勒公式进行展开，但是我们只展开到一阶导数
，把更高阶导数的和看做一个常数量 constant
。
现在我们把对数似然函数按照上式进行展开：
假设 \(\beta^{(t+1)}\) 是对数似然函数的极值点，也就是参数的最优解， \(\beta^t\) 是其附近的一个点。 现在把这个式子进行简单的移项和变换，
显然 \(\ell(\beta^{(t+1)})\) 应该是大于等于 \(\ell(\beta^t)\) 的， 因此有
对上述公式进行移项处理，可得：
我们给参数 \(\beta\) 设置一个初始值，然后通过上式不停的迭代计算新的 \(\beta\) ，\(t\) 表示迭代计算的轮次，直到等号成立的时候，就找到了参数的最优解。
通常我们把一阶导 \(\ell'(\beta^t)\) 称为梯度(gradient)， 公式(8.3.5) 说明只要 \(\beta^{(t+1)}\) 沿着 \(\beta^t\) 的负梯度方向进行移动，我们终将能达到极值点。 注意 \(\frac{constant}{\ell'(\beta^t)}\) 的绝对值的大小影响着前进的速度， 其方向(正负号)决定目标函数是否向着极大值点移动。 所以和下面的公式是等价的，\(\alpha\) 称为学习率(learning rate)，是一个人工设置参数，控制的迭代的速度。
利用 公式(19.2.31) 进行参数迭代求解的方法就称为梯度上升法， 梯度上升法的核心就是让参数变量沿着负梯度的方向前进。 虽然理论上最终一定能到达极值点，但是实际上会受到学习率参数 \(\alpha\) 的影响， 学习率可以理解成每次迭代前进的步长(step size)，步长越大前进的越快，收敛性速度就越快；反之，步长越小，收敛越慢。 但是步长如果大了，就会造成震荡现象，即一步迭代就越过了终点(极值点)，并且在极值点附近往返震荡，永远无法收敛。 为了保证算法能一定收敛，通常会为 \(\alpha\) 设定一个较小的值。 关于 \(\alpha\) 的更多讨论请参考其它资料。
图 8.3.1 梯度下降法中学习率的影响(图片来自网络)¶
8.4. 牛顿法¶
梯度下降法虽然也能收敛到最优解，但是如果学习率设置(通常人工设置)不合理，可能会造成收敛速度太慢或者无法收敛的问题， 其收敛速度难以有效的控制。 现在我们讨论另一中迭代算法，牛顿–拉夫森方法(Newton–Raphson)，一般简称牛顿法。
8.4.1. 算法推导¶
还是从泰勒展开公式开始，让我们考虑二阶泰勒展开：
我们知道目标函数在极值点处的导数应该为 \(0\) ， 所以如果 \(\beta^{(t+1)}\) 是极值点，那么有 \(\ell'(\beta^{(t+1)})=0\) 。我们对 公式(19.2.32) 进行求导，注意 \(\beta^{(t+1)}\) 才是函数未知量， \(\beta_t\) 和 \(\ell(\beta^t)\) 都是已知量。
通过移项可得：
这个迭代等式中，需要同时使用到对数似然函数的一阶导和二阶导数， 二阶偏导数可以在一阶导数的基础上再次求导得到，上一节已经讲过， 对数似然函数的一阶导数又称为得分统计量。
我们对 \(U_j\) 继续求导就是对数似然函数的二阶导数。
对数似然函数的二阶偏导数是一个矩阵，这个矩阵又叫海森矩阵(Hessian matrix) ， 常用符号 \(H\) 表示。牛顿法的迭代公式可以写成如下形式，
和梯度下降法的 公式(19.2.31) 对比下发现，两者非常相似，不同的是牛顿法用 Hessian 矩阵的逆矩阵 \(H(\beta^{(t)})^{-1}\)
替代了学习率参数，避免了需要人工设置学习率的问题。相比梯度下降法，牛顿法收敛速度更快，并且也没有震荡无法收敛的问题。
观察下 公式(8.4.5) ，
GLM 的海森矩阵计算难度是比较大的，为了解决这个问题，
有时候会用海森的矩阵的期望 \(\mathbb{E}[H]\) 替代。
从 公式(8.4.5) 可以看到，海森矩阵是一个关于样本
的函数，所以可以对海森矩阵求关于 \(y\) 的期望。
在参数的迭代过程中使用 \(\mathbb{E}[H]\) 和使用 \(H\) 在参数收敛效果上没有太大区别，二者是类似的，但是 \(\mathbb{E}[H]\) 的计算要简化了很多。 原始海森矩阵 \(H\) 的计算依赖观测样本 \(y_i\) ， 所以通常会把原始海森矩阵称为观测海森矩阵(observed Hessian matrix,OHM) ，他的期望矩阵称为期望海森(expected Hessian matrix,EHM)。
对比下信息矩阵 公式(8.1.15) 和期望海森 公式(8.4.7) ，二者只差一个负号，是相反数的关系，这和我们在 节 2.6 讨论的结论是一致的。
可以看到在 GLM 中，信息矩阵 \(\mathcal{J}\) 可以通过对数似然函数的海森矩阵 \(H\) 得到。
通常把负的 观测 海森矩阵， \(-H\) ， 称为观测信息矩阵(observed information matrix,OIM)，
把负的 期望 海森矩阵， \(- \mathbb{E}[H]\) ， 称为期望信息矩阵(expected information matrix,EIM)。
牛顿法的迭代过程可以用 EIM 代替 OIM 以简化计算过程。
我们这里描述的 Newton–Raphson 算法不支持分散参数 \(\phi\) 的估计， 通常在进行协变量参数 \(\beta\) 的最大似然估计时，认为 \(\phi\) 是已知量。
在 节 2.6 讨论过，参数的最大似然估计估计量是一个统计量，
并且其渐进服从正态分布，其方差可以通过信息矩阵 \(\mathcal{J}\) 计算得到。
最终，Newton–Raphson 提供了如下功能：
为所有
GLM成员模型提供一个参数估计算法。附带产出参数估计量的标准误(standard errors)，可通过信息矩阵得到。
8.4.2. 标准连接函数¶
前文讲过，但模型采用标准连接函数时，得到统计量可以简化。 现在我们看下标准连接函数对牛顿法的影响。 当模型采用标准连接函数时，观测信息矩阵(OIM)会退化成期望信息矩阵(EIM)， 此时在牛顿算法中，两种矩阵是等价的。
根据标准连接函数的定时，当采用标准连接函数时自然参数 \(\theta\) 就等于线性预测器 \(\eta\)，即 \(\theta=\eta\) 此时 \(U\) 可以简化为：
观测海森矩阵是对数似然函数的二阶导数，也是 \(U\) 的一阶导数，因此有
这和 公式(8.4.7) 是一样的。
重要
在 GLM 中，采用标准连接函数时，
观测海森矩阵和期望海森矩阵是相同的，
也就是观测信息矩阵 OIM 和期望信息矩阵 EIM 是相同的。
8.4.3. 迭代初始值的设定¶
要实现 Newton–Raphson 迭代法，
我们必须对参数初始值有一个猜测。
但目前没有用于获得良好参数初值的全局机制，
有一个相对合理的解决方案是，
利用线性预测器中的”常数项系数”获得初始值。
这里的”常数项”指的是线性预测器中截距部分
其中 \(\beta_0\) 就是常数项系数。 如果模型包含常数项，则通常的做法是找到仅包含常数项系数的模型的估计值。 我们令：
然后令对数似然函数的一阶导数 公式(8.1.11) 为 \(0\) ，找到 \(\beta_0\) 的解析解。
通过上式是可以得到 \(\beta_0\) 的一个估计值的。 比如，如果是逻辑回归模型，则有
代入到 公式(8.4.15) 可得：
然后我们就用 \(\beta=(\hat{\beta}_0,0,0,\dots,0)^T\) 作为
Newton–Raphson 算法首次迭代时参数向量的初始值。
如果模型中没有常量项系数，或者我们无法通过解析法求解纯常数项系数模型，则必须使用更复杂的方法，
比如使用搜索方法寻找合理的初始点来开始 Newton-Raphson 算法。
8.5. 迭代重加权最小二乘(IRLS)¶
使用牛顿法对 GLM 中的模型进行参数估计时，
需要把每个模型的对数似然函数通过 \(\beta\) 进行参数化，
然后求出对数似然函数的偏导数，并且在迭代开始前需要给
\(\beta\) 一个初始值，这种方法过于繁琐，
本节我们介绍牛顿法在 GLM 中的一个变种算法，
迭代重加权最小二乘(iteratively reweighted least square,IRLS)算法，
IRLS 算法是``GLM`` 的一个通用型参数估计算法，可用于任意的指数族分布和连接函数，
并且不需要对 \(\beta\) 进行初始化。
8.5.1. 算法推导¶
采用期望海森矩阵的牛顿法的参数迭代等式为
等式两边同时乘以信息矩阵 \(\mathcal{J}\) ，
假设协变量参数 \(\beta\) 的数量是 \(p\) ，则信息矩阵 \(\mathcal{J}\) 是一个 \(p\times p\) 的方阵，其中每个元素 \(\mathcal{J}_{jk}\) 为
仔细观察 \(\mathcal{J}_{jk}\) 的计算公式， 假设有一个 \(N\times N\) 的对角矩阵 ，每个对角元素为
方阵 \(\mathcal{J}\) 就相当于三个矩阵的乘法
这个等式我们先记录下，之后再使用。 现在看下 \(\mathcal{J} \beta\) 的结果是什么。
参数 \(\beta\) 是一个 \(p \times 1\) 的列向量，下标 \(j\) 表示行坐标，下标 \(k\) 表示列坐标。 方阵 \(\mathcal{J}\) 和列向量 \(\beta\) 相乘的计算过程是方阵 \(\mathcal{J}\) 的每个行向量 \(\mathcal{J}_j\) 和列向量 \(\beta\) 进行內积运算，行向量 \(\mathcal{J}_j\) 和 列向量 \(\beta\) 的內积结果为
公式(8.5.2) 的右侧就是两个 \(p\times 1\) 的列向量相加， 每个元素 \(j\) 的计算过程是
我们令
\(Z\) 是一个 \(N \times 1\) 的向量
公式(8.5.2) 的右侧等价于
最终 公式(8.5.2) 等价于
通过移项可以得到参数 \(\beta\) 的迭代公式
其中
\(a(\phi)\) 是分散函数，\(\nu(\mu)\) 是方差函数， \(\frac{\partial \mu}{\partial \eta}\) 是响应函数 \(r\) 的导数， 等价于连接函数 \(g\) 的导数的倒数。 \(\frac{\partial \eta}{\partial \mu}\) 是连接函数 \(g\) 对 \(\mu\) 的导数。 \(W\) 和 \(Z\) 的计算都依赖 \(\eta\) ， 而计算 \(\eta\) 又需要 \(\beta\) ，所以需要迭代的方式更新 \(\beta\) 。
公式(8.5.12) 就是参数向量的更新公式，它在形式上等价于加权的最小二乘法， 其中 \(W\) 相当于权重矩阵，并且每一次迭代都要重新计算 \(W\) ，所以我们把这个算法称为迭代重加权最小二乘法(Iteratively reweighted least square,IRLS)， “reweighted” 指的就是每次迭代重新计算权重矩阵， \(Z\) 被称为工作响应(working response)。
8.5.2. 算法过程¶
收敛性判断
在迭代的过程中，我们可以检查参数 \(\beta\) 的相对变化来决定是否结束算法。
也可以通过相对偏差(deviance)来判断。
关于偏差的概念我们将在下一章详细介绍。
迭代初始值的设定
对比下 Newton–Raphson 算法的参数迭代公式( 公式(8.4.6) ) 和IRLS算法的参数迭代公式( 公式(8.5.12) )， 可以发现IRLS算法并不需要直接在 \(\beta^{(t)}\) 的基础上进行参数迭代， IRLS算法的参数迭代仅仅依赖 \(\mu\) 和 \(\eta\) ，因此与 Newton–Raphson 算法不同的是，IRLS 不需要对参数向量 \(\beta\) 进行初始值的猜测，只需要给 \(\mu\) 和 \(\eta\) 赋予一个初始值即可。
对于二项式分布，可以令 \(\mu_i^{(0)}=k_i(y_i+0.5)/(k_i+1)\) ，\(\eta_i^{(0)}=g(\mu_i^{(0)})\) 。
对于非二项式分布，可以令 \(\mu_i^{(0)}=y_i\) ， \(\eta_i^{(0)}=g(\mu_i^{(0)})\) 。
IRLS算法在更新时，只依赖期望 \(\mu\) 和 线性预测器 \(\eta\) ， 鉴于这一特性，可以使用期望 \(\mu\) 对GLM模型的概率函数进行参数化，而不需要细化到 \(\beta\) ，这可以极大的降低GLM模型概率函数的复杂性。
示例代码
import numpy as np
from typing import Optional, Union
import abc
class Link(abc.ABC):
def link(self, mu: np.ndarray) -> Union[np.ndarray, float]:
"""
连接函数
:param mu:
:return:
"""
pass
def gradient(self, mu: np.ndarray) -> Union[np.ndarray, float]:
"""
连接函数的导数
:param mu:
:return:
"""
pass
def response(self, eta: np.ndarray) -> Union[np.ndarray, float]:
"""
响应函数，连接函数的反函数
:param eta:
:return:
"""
pass
class Distribution(abc.ABC):
name = "default"
def __int__(self, phi=1):
self._phi = phi
def variance(self, mu: np.ndarray) -> Union[np.ndarray, float]:
"""
方差函数
:param mu:
:return:
"""
pass
def phi(self):
"""
分散函数 a(phi)
:return:
"""
return self._phi
class GLM:
def __init__(self, p, link: Link, family: Distribution):
"""
:param p: beta 参数数量
:param link: 连接函数
:param family: 响应变量的分布
:return:
"""
self.beta = np.zeros(shape=(p, 1))
self.link = link
self.family = family
def convergence(self, old_beta, cur_beta):
pass
def init_mu(self, y):
return (y + y.mean()) / 2 # 均值参数初始化
def fit(self, X: np.ndarray, y: np.ndarray):
"""
IRLS算法
:param X:
:param y:
:return:
"""
mu = self.init_mu(y) # 均值参数初始化
eta = self.link.link(mu) # 线性预测器初始化
beta = self.beta.copy()
# 直到收敛
while True:
v = self.family.variance(mu) # 方差函数
g_gradient = self.link.gradient(mu) # 连接函数的导数
W = 1 / (self.family.phi() * v * g_gradient * g_gradient) # 计算权重矩阵
Z = eta + (y - mu) * g_gradient
old_beta = beta
beta = np.invert(X.T * W * X) * X.T * W * Z # 更新参数向量
eta = X * beta # 计算新的线性预测器
mu = self.link.response(eta) # 计算模型预测值/期望参数
if self.convergence(old_beta, beta):
break
self.beta = beta
在 GLM 框架提出之前，各个模型已经被提出并广泛应用了，比如线性回归模型、
逻辑回归模型、泊松回归模型等等，这些模型都是先于 GLM 框架的。
在 GLM 提出前，这些模型都是利用最大似然进行参数估计的，
并且一般都是利用牛顿法进行求解的。
IRLS 算法是伴随着 GLM 诞生的，
要明白的是 IRLS 本身也是建立在最大似然的基础上的。
在 IRLS 提出前，GLM 中的每个模型都需要单独的运用牛顿法求解，
IRLS 是对牛顿法在 GLM 中的一种统一化的抽象。
牛顿法中，每种模型需要单独计算对数似然函数对协变量参数 \(\beta\) 的一、二阶导数，
比较复杂。而 IRLS 建立在 GLM 统一的表达式上，不需要为每个模型单独求参数的导数，
只需要替换相应的连接函数、连接函数导数、方差即可。
重要
IRLS 仍然属于最大似然估计，它是牛顿法在 GLM 中的一种简化。
为了区分，很多资料把采用牛顿法的最大似然估计称为”完全最大似然估计(full maximum likelihood estimation)”。
完全最大似然估计法需要针对每种不同模型单独求解对数似然函数的导数，而 IRLS 不需要。
8.6. 估计量的标准误差¶
在 节 2.7.3 讨论最大似然估计时，讲过参数的最大似然估计量是一个统计量， 而统计量是一个随机量，统计量的概率分布称为抽样分布(sampling distribution)。 期望参数的似然估计量 \(\hat{\mu}\) 的抽样分布是高斯分布。
期望参数的似然估计量 \(\hat{\mu}\) 渐近服从高斯分布，抽样分布的期望值就是参数真实值 \(\mathbb{E}[\hat{\mu}_{ML}] = \mu_{true}\) ，其协方差矩阵是信息矩阵的逆 \(Var(\hat{\mu}_{ML}) = \mathcal{J}^{-1}\) 。 这就意味着均值参数似然估计值的标准误差为
在 GLM 中，期望参数 \(\mu_i\) 和线性预测器 \(\eta_i=\beta^T x_i\) 通过连接函数连接到一起，
协变量参数 \(\beta\) 取代了期望参数 \(\mu_i\) 。
协变量参数 \(\beta\) 的最大似然估计量的抽样分布同样是高斯分布，
这里我们省略证明过程，详细的证明过程可参考 节 10.2 。
在 IRLS 算法的迭代过程中已经计算出了 \(\mathcal{J}(\beta)^{-1}=- \mathbb{E}[H(\beta)]^{-1}=(X^T W X)^{-1}\)
，所以使用 IRLS 算法可以很方便的得到估计量的标准误差。
需要注意的是，IRLS 中使用的是期望信息矩阵（EIM）,
在连接函数是标准连接函数时，它与观测信息矩阵（OIM）相比没有差异，
但当不是标准连接函数时，EIM 的值偏小一些，这会使得 \(SE(\hat{\beta})\) 也偏小，
使我们对参数估计量的标准误差过于乐观，影响我们的判断。
所以当你比较关注参数估计量的标准误差时，建议使用 OIM。
8.7. 分散参数的估计¶
我们已经知道，在所有样本拥有相同分散参数 \(\phi\) 的假设之下， 协变量参数 \(\beta\) 的最大似然估计不会受到 \(\phi\) 的影响。 但这并不意味着分散参数 \(\phi\) 就没有价值了， 首先样本具有相同 \(\phi\) 的假设未必总是成立的， 其次协变量参数 \(\beta\) 的最大似然估计量的标准误差的计算是依赖 \(\phi\) 的。 似然估计量 \(\hat{\beta}_{ML}\) 的标准误是通过信息矩阵 \(\mathcal{J}\) 计算得到的，而 \(\mathcal{J}\) 的计算依赖 \(\phi\) 。
在之前的讨论中，我们都是假设 \(\phi\) 是已知量，通常可以根据人工经验值指定。
然而人工经验不总是靠谱的，很多时候我们需要从实际数据中去探索 \(\phi\) 的合理值。
但是 IRLS 算法并没有提供对 \(\phi\) 的估计，我们需要用一些其它的方法去估计。
最容易想到的方法，就是在得到 \(\beta\) 的最大似然估计值之后，
再次利用最大似然估计对 \(\phi\) 进行估计。
要对 \(\phi\) 进行最大似然估计，就需要对 GLM
的对数似然函数求 \(\phi\) 的导数。
在对数似然函数中， \(c(y_i,\phi)\) 也是包含 \(\phi\) 的，
在 GLM 的不同分布中，它形式是不尽相同的，每种分布模型需要单独去针对 \(\phi\)
求偏导，这种方法比较繁琐，这里我们暂且不表，
本节我们介绍一种简单且常用的方法。
在 GLM 中，估计 \(\phi\) 的最常用方法是利用皮尔逊卡方统计量，
皮尔逊卡方统计量的计算公式为
皮尔逊卡方统计量，顾名思义，它也是一个统计量，并且它的期望值是 \(\mathbb{E}[\chi^2]=N-p-1\) ，\(N\) 是样本的数量，\(p\) 模型输入特征的数量，其实是指协变量参数 \(\beta\) 的数量（不包括截距参数）， \(1\) 表示模型的截距参数。 这里我们假设 \(\chi^2 = N-p-1\) ，\(a(\phi)=\phi\) ，则有
有关皮尔逊卡方统计量的细节在后续章节中会继续讨论，这里可以先记住就可以了。 需要注意的是，利用皮尔逊卡方统计量估计 \(\phi\) 的方法， 同样是建立在所有样本拥有相同 \(\phi\) 的假设之上。