######################################################## 广义线性模型 ######################################################## 我们已经学习讨论了回归模型和分类模型，而分类模型又可以分为生成式模型和判别式模型两大类。 在这些模型中，我们都是在寻找条件概率分布 :math:`p(Y|X)` ， 其中生成式分类模型把输入变量 :math:`X` 和输出变量 :math:`Y` 都看做是随机变量， 然后通过贝叶斯定理的方法求得（后验）条件概率 :math:`p(Y|X)` 。 然而，回归模型(线性回归)和判别式分类模型(比如逻辑回归) 直接为条件概率 :math:`p(Y|X)` 选取一个参数化的概率分布。 本章我们讨论经典线性回归和逻辑回归模型的扩展，在正式讨论前， 我们首先回顾一下线性回顾模型和逻辑回归模型。 **线性回归** 在线性回归模型中，输出变量 :math:`Y` 的值域是实数值，并且假设其服从高斯分布。 此时输出变量 :math:`Y` 的边缘概率分布可以写成： .. math:: p(Y) \sim \mathcal{N}(\mu,\sigma^2) 其中 :math:`\mu` 是均值(期望)参数，表示变量 :math:`Y` 的期望 :math:`\mu=\mathbb{E}[Y]` ， :math:`\sigma^2` 是方差参数，反映变量的波动性。 在回归问题中我们需要对条件概率建模 :math:`p(Y|X)` 进行建模，所谓的条件概率， 就是在变量 :math:`X` 取某个值时(条件下)变量 :math:`Y` 服从什么样的概率分布。 也就是变量 :math:`X` 取不同值时，变量 :math:`Y` 服从"不同的分布"， 这里的"不同分布"通常是不同参数的同种分布。 .. note:: 如果在 :math:`X` 取不同值时，变量 :math:`Y` 服从同一个分布(参数也相同)， 那么说明变量 :math:`X` 和变量 :math:`Y` 是相互独立的，此时有 :math:`p(Y)=P(Y|X)` 。 如果某个输入变量 :math:`X` 与输出变量 :math:`Y` 是相互独立的， 那么这个输入变量就不能作为我们的特征变量，因为它对确定 :math:`Y` 没有任何贡献(影响)。 所以无论是在分类问题还是回归问题中，作为特征的输入变量 :math:`X` 都不能独立于输出变量 :math:`Y` 。 线性回顾模型是通过在边缘概率分布 :math:`p(Y)` 的基础上"引入"变量 :math:`X` 进而得到条件概率分布 :math:`p(Y|X)` 。 我们令变量 :math:`Y` 的期望(均值)值 :math:`\mu` 等于变量 :math:`X` 的线性组合，即： .. math:: \mu = \beta^T x 也就是说在不同的 :math:`x` 值下，变量 :math:`Y` 的均值不同，进而我们就得到了条件概率分布 p(Y|X)： .. math:: p(Y|X) \sim \mathcal{N}(\beta^T x,\sigma^2) 通过参数化 :math:`\mathbb{E}[Y]` 将变量 :math:`X` 引入到边缘概率分布 :math:`p(Y)` 进而得到条件概率分布 :math:`p(Y|X)` 。注意，在经典线性回归模型中，方差 :math:`\sigma^2` 被认为是已知的，并且是和均值独立的， 只有 :math:`\beta` 是需要学习的未知参数。 **逻辑回归** 现在我们回顾一下逻辑回归模型，在逻辑回归模型中，输出变量 :math:`Y` 是伯努利变量，服从伯努利分布： .. math:: p(y) = \mu^y(1-\mu)^{1-y} 在逻辑回归模型中，我们假设均值参数 :math:`\mu` 是一个关于线性组合的 :math:`\beta^Tx` 的sigmoid函数： .. math:: \mu=\mu(x)=\frac{1}{1+e^{-\beta^T x}} 条件概率分布 :math:`p(Y|X)` 就是： .. math:: p(y|x) =\mu(x)^y(1-\mu(x))^{1-y} 我们发现线性回归模型和逻辑回归模型都是通过参数化 :math:`Y` 的期望 :math:`\mu=\mathbb{E}[Y]` ， 进而得到条件概率分布 :math:`p(Y|X)` ，事实上，我们将这种方式进行扩展， 就得到了广义线性模型(generalized linear model,GLM)。 定义 ######################################################## .. _fg_34_1: .. figure:: pictures/34_1.jpg :scale: 70 % :align: center 广义线性模型变量之间的关系 广义线性模型是对经典线性模型的扩展，将输出变量 :math:`Y` 的条件概率分布扩展到指数族分布(指数族的一个子集，不是全部)， :numref:`fg_34_1` 展示了广义线性模型框架下各个变量之间的关系。 - 输入变量 :math:`X` 和系数 :math:`\beta` 组成一个线性关系， :math:`\eta=\beta^T x` ， :math:`\eta` 被称为线性预测器(linear predictor)，:math:`\beta` 是定义的未知参数。 - 通过一个链接函数(link function)将 :math:`\eta` 和变量 :math:`Y` 的均值 :math:`\mu` 关联起来，:math:`\eta=g(\mu)` ， 函数 :math:`g` 称为链接函数。 :math:`g` 的反函数就是激活函数(active function)，:math:`\mu=g^{-1}(\eta)` ，有的资料中也称为响应函数(response function)、均值函数(mean function)。 激活函数可以是线性的，也可以是非线性的， 比如，经典线性回归模型的激活函数为 :math:`\mu=f(\eta)=\eta` ，逻辑回归模型的激活函数为 :math:`\mu=f(\eta)=sigmoid(\eta)` 。 - 在广义线性模型的框架下，响应变量 :math:`Y` 被看做是随机变量，并且其概率分布是指数族分布的一种，:math:`\theta` 是分布的参数。 :math:`\theta` 和 :math:`\mu` 存在一一映射关系，我们用函数 :math:`\psi` 表示这种关系。 显然，一个广义线性模型有三个关键组件： 1. 一个线性预测器 :math:`\eta=\beta^T x` ，被称为系统组件(systematic component)。 2. 一个链接函数 :math:`g` 使得 :math:`\mathbb{E}[Y|X]=\mu=g^{-1}(\eta)` ，链接函数描述系统组件和随机组件之间的关系。 3. 一个指数族分布作为响应变量 :math:`Y` 概率分布 :math:`p(Y;\theta)` ，被称为随机组件(random component)。 指数族分布 =========================================== 我们回顾一下 :numref:`ch_24_1` ，指数族的概率分布的概率密度(质量)函数的一般形式为： .. math:: :label: eq_34_08 p(y|\theta) = \exp \{theta^T T(y) - A(\theta) + S(y)\} 然而在GLM中，通常并不使用上述形式指数族，而是指数族的一个特殊子集，叫做自然指数族(Natural Exponential Family,NEF)， 满足条件 :math:`T(y)=y` 的指数族被称为自然指数族。 .. math:: :label: eq_34_09 p(y|\theta) = \exp \{\theta^T y - A(\theta) + S(y)\} 指数族的这个形式被称为自然形式(natural form)或者规范形式(canonical form)，其中参数 :math:`\theta` 称为自然参数(natural parameter)或者规范参数(canonical parameter)。 自然指数族相对于一般指数族的关键变化就是要求 :math:`T(y)=y` ，因为只有满足这个条件时 :math:`A(\theta)` 的一阶导数才等于 :math:`Y` 的期望 :math:`\mathbb{E}[Y]` 。 在GLM的定义中，我们并不直接使用 :eq:`eq_34_09` 形式，而是在这基础上引入一个额外的参数 :math:`\phi` 。 .. math:: :label: eq_34_EDF p(y|\theta) = \exp \{\frac{\theta y - b(\theta)}{a(\phi)} + c(y,\phi)\} 这种形式的指数族通常被称为指数分散族(exponential dispersion family,EDF)， 参数 :math:`\phi` 称为分散参数(dispersion parameter)， :math:`a(\phi)` 称为分散函数(dispersion function)， 分散参数 :math:`\phi` 和分布的方差有关。 并不是所有的指数族分布都存在 :math:`\phi` ，比如对于泊松分布、二项式分布来说，:math:`\phi=1` 。 在EDF中，:math:`\theta` 不再是向量，而是一个标量，仍然称为 canonical parameter， **参数** :math:`\theta` **是和** :math:`Y` **的均值** :math:`\mu` **相关的，** :math:`\theta` 和 :math:`\mu` 之间存在一一映射的函数关系 :math:`\theta=\psi(\mu)` 。 换句话说，:math:`\theta` 和 :math:`\mu` 可以互相转化。 通常在GLM中，只有参数 :math:`\theta` 作为模型的未知参数，此时称为单参数模。 单参数模型指的是模型中只有 :math:`\theta` 是未知参数，而 :math:`\phi` 是已知的， **分散参数(dispersion parameter)** 在最初的GLM论文中(Nelder and Wedderburn, 1972)把 :math:`a(\phi)` 称为比例因子(scale factor)， 并且没有给参数 :math:`\phi` 单独命名。后来在1974年 Royal Statistical Society 发布了首个 GLM的软件工具包(Generalized Linear Interactive Modelling,GLIM)， 在GLIM中把 :math:`a(\phi)` 定义成： .. math:: a(\phi) = \frac{\phi}{w} 其中 :math:`w` 是先验权重(prior weight)，并且把 :math:`\phi` 称为尺度参数(scale parameter)， 这就是导致了对 :math:`\phi` 命名产生了歧义。 因为"scale"这个词在统计学还有其它用法，容易产生混淆， 所以在1980s(McCullagh and Nelder)初版的GLM书籍中， 把 :math:`\phi` 命名为"dispersion parameter"， 后来也就沿用了这种叫法。 但是由于GLIM流行了很久，导致"scale"的叫法还存在很多资料中。 在很多GLM的工具包中，都会把 :math:`a(\phi)` 定义成如下形式： .. math:: a(\phi)_n = \frac{\phi}{w_n} 其中 :math:`w_n` 是观测样本的权重，不同的样本可以拥有不同的权重值， 比如在利用ML进行参数估计时，对于某些样本设置成 :math:`w_n=0` ，这就相当于抛弃了这些样本。 **累积函数(cumulant function)** 我们知道在 :eq:`eq_34_09` 的指数族形式中 :math:`A(\theta)` 称为累积函数(cumulant function)， 可以用 :math:`A(\theta)` 的导数求出分布的矩，一阶导数是分布的期望， 二阶导数是分布的方差。 然而在GLM中我们使用的是 :eq:`eq_34_EDF` 的形式， 其中 :math:`b(\theta)` 本质上就是 :math:`A(\theta)` ，二者关系是： .. math:: A(\theta) = \frac{ b(\theta) }{a(\phi)} 所以我们同样把 :math:`b(\theta)` 被称为累积函数(cumulant function)， 并且它同样和分布的矩(moments)有关。 .. math:: \mathbb{E}[Y] = b'(\theta)=\mu .. math:: :label: eq_34_20 Var(Y) = A''(\theta)=a(\phi)b''(\theta) **方差结构** 在EDF(指数分散族，Exponential Dispersion Family)中， 分布的方差可以表示成两部分的乘积（ :eq:`eq_34_20` ）， 一部分是分散函数 :math:`a(\phi)` ，另一部分是累计函数的二阶导数 :math:`b''(\theta)` 实际上，:math:`b''(\theta)` 是一个关于 :math:`\mu` 函数， 首先 :math:`b(\theta)` 是一个关于 :math:`\theta` 的函数， 其二阶导数自然也是一个关于 :math:`\theta` 的函数(或者是一个常数)。 然而自然参数 :math:`\theta` 和均值参数 :math:`\mu` 存在一一对应关系，所以一定可以把 :math:`\theta` 替换成 :math:`\mu` ，所以 :math:`b''(\theta)` 一定可以写成一个关于 :math:`\mu` 的函数， 因此我们定义 :math:`\nu(\mu)=b''(\theta)` 。 .. math:: Var(Y) =\sigma^2 = b''(\theta) a(\phi) = \nu(\mu)a(\phi) 函数 :math:`\nu(\mu)=b''(\theta)` 称为方差函数(variance function)，其将分布的均值参数 :math:`\mu` 和分布的方差关联在一起。 如果其值一个常数值，说明均值和方差是独立无关的；反之，如果其最终是 :math:`\mu` 的函数，说明均值和方差是相关联的。 在统计中，方差函数是一个平滑函数，它把随机量的方差描述成其均值的函数。 在高斯分布中， :math:`b''(\theta)=1` ，所以方差和均值是相互独立的， 对于其他分布，这是不成立的，这使得高斯分布是特例。 分散参数 :math:`\phi` 也影响着分布的方差， 参数 :math:`\theta` 和 :math:`\phi` 本质上是位置(locate)和比例(scale)参数， 位置参数反映数据的均值，比例参数反映数据方差。 当然对于很多分布， :math:`a(\phi)=1` 。 在经典线性回归模型中，输入特征数据 :math:`x` 通过线性组合 :math:`\eta=\beta^T x` 影响着响应变量 :math:`Y` (高斯分布) 的均值 :math:`\mu=\eta=\beta^T x` ， 所有的观测样本共用参数 :math:`\beta` (对于任意 :math:`x` ，都是同样的 :math:`\beta` 值)， 当 :math:`x` 不同时， 高斯变量 :math:`Y` 拥有不同的均值 :math:`\mu` ， 通过这种方式实现了条件概率分布 :math:`p(Y|X)` 的表达。 但是对于高斯变量 :math:`Y` 的方差参数 :math:`\sigma^2` 并没有假设为未知参数，而是假设其为已知的值， 并且对于任意的观测样本 :math:`x` 都是一样的值。 然而，在GLM的框架下，是可以允许不同观测样本有不同的方差， 而这是通过 :math:`a(\phi)` 实现的。 此时函数 :math:`a(\phi)` 通常被定义成如下的形式： .. math:: a(\phi) = \frac{\phi}{w_n} 通常对于所有的观测样本来说，:math:`\phi` 是一个相同的，而 :math:`w` 可以根据不同的观测样本取不同的值， 下标 :math:`n` 表示样本编号。 :math:`w` 被称为先验权重(prior weight)，通常是根据额外的先验信息确定的。 如果所有观测样本具有相同的方差假设，那么 :math:`w` 值通常就是1；反之，:math:`w` 可以是和样本相关的， 不同的样本采用不同的值。 .. csv-table:: 常见分布的方差函数 :header: "分布","方差函数", "约束", "导数 :math:`\\partial \\nu(\\mu) / \\partial\\mu`" "Gaussian", ":math:`1`", ":math:`\left\{ \begin{array}{lr} \mu \in \mathcal{R} \\ y \in \mathcal{R} \end{array} \right .`", ":math:`0`" "Bernoulli", ":math:`\mu(1-\mu)`", ":math:`\left\{ \begin{array}{lr} 0<\mu<1 \\ 0 \le y \le 1 \end{array} \right .`", ":math:`1-2\mu`" "Binomial(k)", ":math:`\mu(1-\mu/k)`", ":math:`\left\{ \begin{array}{lr} 0<\mu0 \\ y \ge 0 \end{array} \right .`", ":math:`1`" "Gamma", ":math:`\mu^2`", ":math:`\left\{ \begin{array}{lr} \mu >0 \\ y > 0 \end{array} \right .`", ":math:`2\mu`" "Inverse Gaussian", ":math:`\mu^3`", ":math:`\left\{ \begin{array}{lr} \mu >0 \\ y > 0 \end{array} \right .`", ":math:`3\mu^2`" "Negative binomial(:math:`\alpha`)", ":math:`\mu+\alpha\mu^3`", ":math:`\left\{ \begin{array}{lr} \mu >0 \\ y \ge 0 \end{array} \right .`", ":math:`1+2\alpha\mu`" "Power(:math:`k`)", ":math:`\mu^k`", ":math:`\left\{ \begin{array}{lr} \mu >0 \\ k \ne 0,1,2 \end{array} \right .`", ":math:`k\mu^{k-1}`" "Quasi", ":math:`\nu(\mu)`", "", ":math:`\frac{\partial \nu(\mu)}{\partial \mu}`" 链接函数 ======================================= 链接函数 :math:`g` 是用来链接线性预测器 :math:`\eta` 和均值 :math:`\mu` 的，要求 **必须是连续可微，并且可逆的。** 原则上，任何单调的连续可微函数都可以，但是对于标准GLM，有一些方便且通用的选择。 在高斯线性模型中，链接函数是恒等函数 :math:`\eta=g(\mu)=\mu` 。 在泊松分布中，均值 :math:`\mu` 必须是正的，所以 :math:`\eta=\mu` 不再适用，因为 :math:`\eta=\beta^Tx` 的取值范围值整个实数域。对于泊松分布，标准的链接函数选择是对数函数 :math:`\eta=log \mu` ，此时 :math:`\mu=e^{\eta}` 确保了 :math:`\mu` 为正数。 **链接函数本质上，就是把实数域范围的** :math:`\eta` **转换到特定分布合法的** :math:`\mu` **值空间上。** .. topic:: 规范链接(canonical link) 当链接函数使得 :math:`\eta=\theta` 时，称为规范链接(canonical link)函数。 实际上规范链接函数满足 :math:`\eta=g(\mu)=\psi(\mu)=\theta` ， 换句话说，对于一个特定的指数族分布，其规范链接函数为 :math:`g=\psi` 。 使用规范链接函数可以带来很多统计属性，最直接的就是可以简化参数估计的计算过程。 .. csv-table:: 常见链接函数 :header: "名称","链接函数", "激活函数(反链接)", " :math:`\\mu` 的空间" "Identity", ":math:`\eta=\mu`", ":math:`\mu=\eta`", ":math:`\mu \in \mathcal{R}`" "Logit", ":math:`\eta=\ln\{\mu/(1-\mu)\}`", ":math:`\mu=e^\eta/(1+e^\eta)`", ":math:`\mu \in (0,1)`" "Log", ":math:`\eta=\ln(\mu)`", ":math:`\mu=e^\eta`", ":math:`\mu >0`" "Negative binomial( :math:`\alpha`)", ":math:`\eta=\ln\{\mu/(\mu+1/\alpha)\}`", ":math:`\mu=e^\eta/\{ \alpha(1-e^\eta)\}`", ":math:`\mu >0`" "Log-complement", ":math:`\eta=\ln(1-\mu)`", ":math:`\mu=1-e^\eta`", ":math:`\mu <1`" "Log-log", ":math:`\eta=-ln \{- \ln(\mu)\}`", ":math:`\mu=\exp\{-\exp(-\eta)\}`", ":math:`\mu \in (0,1)`" "Complementary log-log", ":math:`\eta=ln \{- \ln(1-\mu)\}`", ":math:`\mu=1-\exp\{-\exp(\eta)\}`", ":math:`\mu \in (0,1)`" "Probit", ":math:`\eta=\Phi^{-1}(\mu)`", ":math:`\mu=\Phi(\eta)`", ":math:`\mu \in (0,1)`" "Reciprocal", ":math:`\eta=1/\mu`", ":math:`\mu=1/\eta`", ":math:`\mu \in \mathcal{R}`" "Power(:math:`\alpha=-2`)", ":math:`\eta=1/\mu^2`", ":math:`\mu=1/\sqrt{\eta}`", ":math:`\mu >0`" "Power(:math:`\alpha`) :math:`\left\{ \begin{array}{lr}\alpha \ne 0\\ \alpha=0 \end{array} \right .`", ":math:`\eta=\left\{ \begin{array}{lr}\mu^\alpha \\ \ln(\mu) \end{array} \right .`", ":math:`\mu=\left\{ \begin{array}{lr}\eta^{1/\alpha} \\ \exp(\eta) \end{array} \right .`", ":math:`\mu \in \mathcal{R}`" "Odds power(:math:`\alpha`) :math:`\left\{ \begin{array}{lr}\alpha \ne 0\\ \alpha=0 \end{array} \right .`", ":math:`\eta=\left\{ \begin{array}{lr} \frac{\mu/(1-\mu)^\alpha-1}{\alpha} \\ \ln \left( \frac{\mu}{1-\mu} \right) \end{array} \right .`", ":math:`\mu=\left\{ \begin{array}{lr} \frac{(1+\alpha\eta)^{1/\alpha}}{1+(1+\alpha\eta)^{1/\alpha}} \\ \frac{e^\eta}{1+e^\eta} \end{array} \right .`", ":math:`\mu \in (0,1)`" .. csv-table:: 链接函数的导数 :header: "名称","链接函数", "一阶导数 :math:`\\triangle=\\partial \\eta /\\partial \\mu`", "二阶导数" "Identity", ":math:`\eta=\mu`", ":math:`1`", ":math:`0`" "Logit", ":math:`\eta=\ln\{\mu/(1-\mu)\}`", ":math:`1/\{\mu(1-\mu)\}`", ":math:`(2\mu-1)\triangle^2`" "Log", ":math:`\eta=\ln(\mu)`", ":math:`1/\mu`", ":math:`-\triangle^2`" "Negative binomial( :math:`\alpha`)", ":math:`\eta=\ln\{\alpha\mu/(1+\alpha\mu)\}`", ":math:`1/(\mu+\alpha\mu^2)`", ":math:`-\triangle^2(1+2\alpha\mu)`" "Log-complement", ":math:`\eta=\ln(1-\mu)`", ":math:`-1/(1-\mu)`", ":math:`-\triangle^2`" "Log-log", ":math:`\eta=-ln \{- \ln(\mu)\}`", ":math:`-1/\{\mu\ln(\mu)\}`", ":math:`\{1+\ln(\mu)\}\triangle^2`" "Complementary log-log", ":math:`\eta=ln \{- \ln(1-\mu)\}`", ":math:`\{(\mu-1)\ln(1-\mu)\}^{-1}`", ":math:`-\{1+\ln(1-\mu)\}\triangle^2`" "Probit", ":math:`\eta=\Phi^{-1}(\mu)`", ":math:`1/\phi\{\Phi^{-1}(\mu)\}`", ":math:`\eta\triangle^2`" "Reciprocal", ":math:`\eta=1/\mu`", ":math:`-1/\mu^2`", ":math:`-2\triangle / \mu`" "Power(:math:`\alpha=-2`)", ":math:`\eta=1/\mu^2`", ":math:`-2/\mu^3`", ":math:`-3\triangle / \mu`" "Power(:math:`\alpha`) :math:`\left\{ \begin{array}{lr}\alpha \ne 0\\ \alpha=0 \end{array} \right .`", ":math:`\eta=\left\{ \begin{array}{lr}\mu^\alpha \\ \ln(\mu) \end{array} \right .`", ":math:`\left\{ \begin{array}{lr} \alpha \mu ^{\alpha-1} \\ 1/\mu \end{array} \right .`", ":math:`\left\{ \begin{array}{lr} (\alpha-1)\triangle/\alpha \\ -\triangle^2 \end{array} \right .`" "Odds power(:math:`\alpha`) :math:`\left\{ \begin{array}{lr}\alpha \ne 0\\ \alpha=0 \end{array} \right .`", ":math:`\eta=\left\{ \begin{array}{lr} \frac{\mu/(1-\mu)^\alpha-1}{\alpha} \\ \ln \left( \frac{\mu}{1-\mu} \right) \end{array} \right .`", ":math:`\left\{ \begin{array}{lr} \frac{\mu^{\alpha-1}}{(1-\mu)^{\alpha+1}} \\ \frac{1}{\mu(1-\mu)} \end{array} \right .`", ":math:`\left\{ \begin{array}{lr} \triangle\left(\frac{1-1/\alpha}{1-\mu} +\alpha+1\right) \\ \mu\triangle^2 \end{array} \right .`" .. csv-table:: 激活函数的导数 :header: "链接函数名称", "激活函数(反链接)","一阶导数 :math:`\\triangle=\\partial \\mu /\\partial \\eta`", "二阶导数" "Identity", ":math:`\mu=\eta`", ":math:`1`", ":math:`0`" "Logit", ":math:`\mu=e^\eta/(1+e^\eta)`", ":math:`\mu(1-\mu)`", ":math:`\triangle(1-2\mu)`" "Log", ":math:`\mu=e^\eta`", ":math:`\mu`", ":math:`\triangle`" "Negative binomial( :math:`\alpha`)", ":math:`\mu=e^\eta/\{ \alpha(1-e^\eta)\}`", ":math:`\mu+\alpha\mu^2`", ":math:`\triangle(1+2\alpha\mu)`" "Log-complement", ":math:`\mu=1-e^\eta`", ":math:`\mu-1`", ":math:`\triangle`" "Log-log", ":math:`\mu=\exp\{-\exp(-\eta)\}`", ":math:`-\mu\ln(\mu)`", ":math:`\triangle\{ 1+\ln(\mu)\}`" "Complementary log-log", ":math:`\mu=1-\exp\{-\exp(\eta)\}`", ":math:`(\mu-1)\ln(1-\mu)`", ":math:`\triangle\{1+\ln(1-\mu)\}`" "Probit", ":math:`\mu=\Phi(\eta)`", ":math:`\phi(\eta)`", ":math:`-\triangle \eta`" "Reciprocal", ":math:`\mu=1/\eta`", ":math:`-\mu^2`", ":math:`-2\triangle\mu`" "Power(:math:`\alpha=-2`)", ":math:`\mu=1/\sqrt{\eta}`", ":math:`-\mu^3/2`", ":math:`3\triangle^2 / \mu`" "Power(:math:`\alpha`) :math:`\left\{ \begin{array}{lr}\alpha \ne 0\\ \alpha=0 \end{array} \right .`", ":math:`\mu=\left\{ \begin{array}{lr}\eta^{1/\alpha} \\ \exp(\eta) \end{array} \right .`", ":math:`\left\{ \begin{array}{lr} \frac{1}{\alpha} \mu^{1-\alpha} \\ \mu \end{array} \right .`", ":math:`\left\{ \begin{array}{lr} \triangle(1/\alpha -1)/\mu^\alpha \\ \triangle \end{array} \right .`" "Odds power(:math:`\alpha`) :math:`\left\{ \begin{array}{lr}\alpha \ne 0\\ \alpha=0 \end{array} \right .`", ":math:`\mu=\left\{ \begin{array}{lr} \frac{(1+\alpha\eta)^{1/\alpha}}{1+(1+\alpha\eta)^{1/\alpha}} \\ \frac{e^\eta}{1+e^\eta} \end{array} \right .`", ":math:`\left\{ \begin{array}{lr} \frac{\mu(1-\mu)}{1+\alpha\eta} \\ \mu(1-\mu) \end{array} \right .`", ":math:`\left\{ \begin{array}{lr} \triangle \left ( 1-2\mu-\frac{\alpha}{1+\alpha\eta} \right) \\ \triangle(1-2\mu) \end{array} \right .`" 当无法合理地假设数据是正态分布的或者响应变量的结果集有限集时，传统的线性回归模型是不合适的。 此外，在许多情况下，传统线性回归模型的同方差假设是站不住脚的，此时传统线性回归模型也是不合适的。 GLM允许对传统线性回归模型进行扩展，以突破这些限制。 GLM框架通过链接函数把线性预测变量 :math:`\eta=\beta^Tx` 映射到服从指数族分布的响应变量的均值参数 :math:`\mu` 。并且我们能够开发一种适用于所以GLM框架下模型的参数估计算法， 所以我们可以支持对诸如logit，probit和Poisson之类的有用统计模型的参数估计。 例子 ============================================ **高斯分布** 高斯分布的概率密度函数为： .. math:: f(y) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \{ -\frac{1}{2}\frac{(y-\mu)^2}{\sigma^2} \} 把其改写成指数分散族的形式： .. math:: f(y) = \exp \{ \frac{y\mu-\frac{1}{2}\mu^2}{\sigma^2} - \frac{y^2}{2\sigma^2} - \frac{1}{2} \ln (2\pi\sigma^2) \} 和 :math:`eq_34_EDF` 进行对比，各个标准项为： .. math:: \theta &=\mu b(\theta) &= \frac{1}{2}\mu^2 a(\phi) &= \sigma^2 高斯分布的均值和方差为： .. math:: \mathbb{E}[Y]=b'(\theta) = \mu Var(Y) = b''(\theta)a(\phi) = \sigma^2 对于高斯分布来说，方差和均值是独立无关的。 **伯努利分布** .. csv-table:: 常见GLM表(1) :header: "","Normal(Gaussian) :math:`N(\\mu,\\sigma^2)`", "Bernoulli :math:`B(\\mu)`", "Binomial :math:`B(N,\\mu)`" "Range of y", "real: :math:`(-\infty,+\infty)`", ":math:`\{0,1\}`", ":math:`\{0,\dots,N\}`" "f(y)", :math:`\frac{1}{\sqrt{2\pi\sigma^2}}\exp \{ -\frac{(y-\mu)^2}{2\sigma^2} \}`, :math:`\mu^y(1-\mu)^{1-y}`,:math:`\binom{N}{y}\mu^y(1-\mu)^{N-y}` "EDF",:math:`\exp\{\frac{\mu y-\frac{\mu^2}{2}}{\sigma^2}-\frac{y^2}{2\sigma^2}-\frac{\ln 2\pi\sigma^2}{2}\}`, :math:`\exp\{y \ln \frac{\mu}{1-\mu} + \ln(1-\mu)\}`,:math:`\exp\bigg[ \frac{y \ln(\frac{\mu}{1-\mu}) + \ln(1-\mu)}{1/N} + \ln({N \choose y})\bigg]` ":math:`\theta=\psi(\mu)`", :math:`\theta=\mu`, :math:`\theta=\ln \left ( \frac{\mu}{1-\mu} \right )=logit(\mu)`,:math:`\theta=\ln \left ( \frac{\mu}{1-\mu} \right )` ":math:`\mu=\psi^{-1}(\theta)`", :math:`\mu=\theta`, :math:`\mu=\frac{1}{1+e^{-\theta}}=sigmoid(\theta)`,:math:`\mu=\frac{1}{1+e^{-\theta}}` ":math:`b(\theta)`",:math:`\frac{\theta^2}{2}`,:math:`\ln(1+e^{\theta})`,:math:`\ln(1+e^{\theta})` ":math:`b(\mu)`",:math:`\frac{\mu^2}{2}`,:math:`-\ln(1-\mu)`,:math:`-\ln(1-\mu)` "Link name", Identity, Logit,Logit "Link function", :math:`\eta=\mu`, :math:`\eta=\ln \left( \frac{\mu}{1-\mu} \right)`,:math:`\eta=\ln \left( \frac{\mu}{1-\mu} \right)` "Mean function", :math:`\mu=\eta`, :math:`\mu=\frac{1}{1+e^{-\eta}}`,:math:`\mu=\frac{N}{1+e^{-\eta}}` ":math:`\nu(\mu)=b''(\theta)`", 1, :math:`\mu(1-\mu)`,:math:`\mu(1-\mu)` ":math:`a(\phi)`", :math:`\sigma^2`,1,:math:`\frac{1}{N}` .. csv-table:: 常见GLM表(2) :header: "", "Categorical :math:`Cat(K,\\mu)`","Poisson :math:`Poisson(\\mu)`" "Range of y", ":math:`\{1,\dots,K\}`","integer :math:`0,1,2,\dots`" "f(y)", :math:`\prod_{k}\mu_k^{y_k}`, :math:`\exp\{y\ln \mu - \ln\mu\}` "EDF",:math:`\exp \left \{ \sum_{k=1}^{K-1} x_k \ln \left ( \frac{\mu_k}{ \mu_K} \right )+ \ln \left (1-\sum_{k=1}^{K-1} \mu_k \right ) \right \}`,:math:`\exp\{y\ln \mu - \ln\mu\}` ":math:`\theta=\psi(\mu)`", :math:`\theta_k=\ln \left ( \frac{\mu_k}{\mu_K} \right )`,:math:`\theta=\ln \mu` ":math:`\mu=\psi^{-1}(\theta)`",:math:`\mu_k = \frac{e^{\theta_k}}{\sum_{j=1}^K e^{\theta_j}}`,:math:`\mu=e^{\theta}` ":math:`b(\theta)`",:math:`\ln \left ( \sum_{k=1}^K e^{\theta_k} \right )`,:math:`e^{\theta}` ":math:`b(\mu)`",:math:`- \ln \left (1-\sum_{k=1}^{K-1} \mu_k \right )`,:math:`\ln\mu` "Link name", Logit "Link function", :math:`\eta_k=\ln \left( \frac{\mu_k}{\mu_K} \right)` "Mean function", :math:`\mu_k=\frac{e^{\eta_k}}{\sum_k e^{\eta_k}}` ":math:`\nu(\mu)=b''(\theta)`", :math:`\mu_k(1-\mu_k)`,\mu ":math:`a(\phi)`", 1,1 参数估计 ######################################################## 本节我们介绍两种GLM模型的参数估计算法， 我们将统一的以指数族的形式展现算法过程，这样适用于指数族中的所有具体分布， 我们的目标是让大家对GLM的基础理论有个全面的了解，同时我们会着重强调算法成立的假设及其一些限制。 在本章的后几节我们将说明该理论在特定指数族成员上的运用。 传统上，对于单参数的指数族分布可以运用梯度下降法和牛顿法进行参数估计， 梯度下降法的优点是算法实现简单，缺点是收敛速度不如牛顿法。 梯度下降法和牛顿法在形式上是非常相似的，二者都是沿着目标函数的负梯度方向寻找最优解， 不同的是传统梯度下降法利用一阶导数，而牛顿法利用二阶导数，牛顿法相对于梯度下降法收敛速度会更快， 但是由于二阶导数的引入也使得牛顿法的计算复杂度增加很对，甚至很多时候无法计算。 所以这里我们同时会介绍牛顿法的一个变种算法，迭代重加权最小平方法(iteratively reweighted least squares,IRLS)。 GLM框架下的模型都可以以统计的形式运用IRLS算法进行参数估计，这是GLM非常有吸引力的一点。 在之后的参数估计算法讨论中，我们都是假设在一个满足独立同分布(IID)的观测样本集上进行参数估计， 并且这个样本集是指数族分布的观测样本集。所有样本的联合概率通常可以被称为样本集的似然函数， 独立同分布的样本集的联合概率函数可以写成： .. math:: :label: eq_34_21 f(y;\theta,\phi)=\prod_{n=1}^N f(y_n;\theta,\phi) 似然函数可以理解为： .. math:: :label: eq_34_22 L(\theta,\phi;y)=\prod_{n=1}^N f(\theta,\phi;y_n) :eq:`eq_34_21` 和 :eq:`eq_34_22` 的区别在于， 前者是一个概率密度(质量)函数，其是在给定 :math:`\theta,\phi` 的条件下关于未知量 :math:`y` 的函数； 后者是一个似然函数，其表达是在给定数据 :math:`y` 的条件下关于未知参数 :math:`\theta,\phi` 的函数。 因为联合似然函数是一个函数的连乘形式，所以通常我们会对其进行一个对数转换(log-transform)进而得到一个连加的形式， 连加的形式更方便进行计算。加了对数的似然函数被称为对数似然函数，:math:`\ell` ，这个函数也是最大似然估计(ML)的核心，我们从指数族的对数似然函数开始。 .. note:: 连乘变成连加有两个好处，(1)更容易求导和极大化操作；(2)似然函数是概率连乘，而概率都是小于1的，大量小于1的数字连乘产生更小的数字， 甚至趋近于0，而计算机的浮点数精度是通常无法处理这么小的数字的，所以加对数更方便计算机进行数值处理。 .. math:: \ell=\sum_{n=1}^N \left \{ \frac{y_n\theta_n - b(\theta_n)}{a(\phi)} + c(y_i,\phi) \right \} :math:`\theta` 是规范参数(canonical parameter )； :math:`b(\theta)` 是累积量，它描述的是分布的矩(moment)； :math:`\phi` 是分散参数(dispersion parameter)，它是比例或辅助参数； :math:`c(\cdot)` 是归一化项。 归一化项不是 :math:`\theta` 的函数，而是简单地缩放基础密度函数的范围使得整个函数的积分（或求和）为1。 在上一节我们讨论过，指数族分布的期望与方差可以通过 :math:`b(\theta)` 的导数求得。 .. math:: \mu &= \mathbb{E}[Y] = b'(\theta) Var(Y) &= b''(\theta) a(\phi) 并且我们知道，均值参数 :math:`\mu` 和规范参数 :math:`\theta` 是存在一个可逆的函数关系的，也就是说 :math:`\mu` 可以看做是关于 :math:`\theta` 的一个函数，反之，:math:`\theta` 也可看做是一个关于 :math:`\mu` 的函数。 基于这个事实，我们可以把 :math:`b''(\theta)` 看做是一个关于 :math:`\mu` 的函数，记作 :math:`\nu(\mu)` 。 .. math:: b''(\theta) \triangleq \nu(\mu) 因此，方差 :math:`Var(Y)` 就可以被看成是函数 :math:`\nu(\mu)` 和分散函数 :math:`a(\phi)` 的乘积，通常我们把 :math:`\nu(\mu)=b''(\theta)` 称为方差函数(variance function)， **注意：虽然叫方差函数，但方差函数的值不是方差本身。** 有时 :math:`b''(\theta)` 会是一个常数量(constant)，比如高斯分布，此时分布的方差为： .. math:: Var(Y)=constant \times a(\phi) 这时，分布的方差就不会受到均值的影响了。 另外方差函数 :math:`\nu(\mu)` 可以通过简单方式求得。 .. math:: \nu(\mu) = b''(\theta) =(b'(\theta))'= (\mu(\theta))' = \frac{\partial \mu}{\partial \theta} 显然，当 :math:`\mu` 与 :math:`\theta` 之间的映射函数是线性函数时，一阶偏导 :math:`\frac{\partial \mu}{\partial \theta}` 就是一个常数值。另外，我们知道反函数的导数就等于原函数导数的倒数，所以有： .. math:: \frac{\partial \theta}{\partial \mu} = \frac{1}{\nu(\mu)} 在GLM框架下，输入变量 :math:`X` 和其系数 :math:`\beta` 组成一个线性预测器 :math:`\eta=\beta^Tx` 。 :math:`\eta` 和分布的均值(期望)通过链接函数(已知的)连接在一起 。 .. math:: \eta &=\beta^Tx = g(\mu) \mu &= g^{-1}(\eta) 其中 :math:`\beta` 和 :math:`x` 都是一个向量，:math:`\beta^Tx = \sum_j \beta_j x_j` ， 因此有： .. math:: \frac{\partial \eta}{\partial \beta_j} = x_j 线性预测器 :math:`\eta` 的值空间并没有特别的限定， 其值空间可以是整个实数域 :math:`\eta \in R` 。 **因此，链接函数的一个目的就是将线性预测器的值映射到输出变量的范围。** 映射到输出变量的范围可以确保和指数族特定成员分布相关。 一个特殊的链接函数是恒等函数， :math:`\eta=\theta` ， 恒等函数的链接函数被称为规范链接(canonical link)或者自然链接(natural link)。 现在让我们回到似然估计，似然估计的一个关键过程就是要对数似然函数的求导。 经过以上的铺垫，我们利用求导的链式法则对GLM模型的对数似然函数进行求导。 注意，我们的参数 :math:`\beta` 是一个向量参数，我们只需要对其中的一项进行求导即可。 .. math:: :label: eq_34_jac \frac{ \partial \ell}{\beta_j} &= \sum_{n=1}^N \left ( \frac{\partial \ell_n}{\partial \theta_n} \right ) \left ( \frac{\partial \theta_n}{\partial \mu_n} \right ) \left ( \frac{\partial \mu_n}{\partial \eta_n} \right ) \left ( \frac{\partial \eta_n}{\partial \beta_j} \right ) &= \sum_{n=1}^N \left \{ \frac{y_n-b'(\beta_n)}{a(\phi)} \right \} \left \{ \frac{1}{\nu(\mu_n)} \right \} \left ( \frac{\partial \mu}{\partial \eta} \right )_n x_{jn} &= \sum_{n=1}^N \frac{y_n-\mu_n}{a(\phi) \nu(\mu_n) } \left ( \frac{\partial \mu}{\partial \eta} \right )_n x_{jn} 其中 :math:`n` 是观测样本的编号，:math:`x_{jn}` 表示第 :math:`n` 条观测样本的第 :math:`j` 列特征值，其值已知。 :math:`y_n` 是当前观测样本的标签值，也是已知的。 :math:`\mu_n=g^{-1}(\eta_n)=g^{-1}(\beta^Tx_n)` 是链接函数的反函数的值，代入样本值 :math:`x_n` 和当前的参数值 :math:`\beta` 后可以直接算出，方差函数 :math:`\nu(\mu_n)` 是关于 :math:`\mu_n` 的函数，因此也可以算出。 :math:`a(\phi)` 是已知的， :math:`\frac{\partial \mu}{\partial \eta} ` 是函数 :math:`g^{-1}(\eta_n)` 关于 :math:`\eta_n` 的导数， 在确定了链接函数的形式后也是可以算出的。 :eq:`eq_34_jac` 是GLM标准指数族形式下对数似然函数的一阶偏导数，GLM框架下的任意模型都可以按照此公式计算偏导数， 只需要按照特定的分布和链接函数替换相应组件即可。 回顾一下线性回归的章节，要想求得参数的估计值，我们可以令对数似然函数的一阶导数为0，得到正规方程(normal equations)， 然后解正规方程的方式得到解析解。 然而，正规方程并不是一定能求解的，需要满足一些限制条件才行， 比如链接函数必须是规范链接函数(canonical link)以及满足矩阵可逆(如果不是很理解可以回顾一下线性回归的章节- :numref:`ch_29` )。 所以解析解的方式并不具备通用性，我们需要采用更一般的方法，迭代法。 迭代法又可以简单分为一阶导(梯度下降法系列)和二阶导(牛顿法系列)，实际这两种都可以通过泰勒展开(Taylor explanation)公式进行证明。 这里我们简单回顾一下泰勒展开定理。 注意，本书讨论的迭代求解算法默认目标函数都是凸函数，也就是函数有唯一的极值点。 关于非凸函数以及带约束的优化问题，请读者参考其它资料。 .. topic:: 泰勒展开 设 :math:`n` 是一个正整数。如果定义在一个包含 :math:`x_0` 的区间上的函数f在 :math:`x_0` 处 :math:`n+1` 次可导，那么对于这个区间上的任意x，都有 .. math:: f(x)_{Taylor} &= \sum_{n=0}^{\infty} \frac{f^{(n)}(x_0)}{n!} \times (x - x_0)^n &= f(x_0) + \frac{f'(x_0)}{1!}(x-x_0) + \frac{f^{(2)}(x_0)}{2!}(x-x_0)^2+ \cdots + \frac{f^{(n)}(x_0)}{n!}(x-x_0)^n + R_n(x) .. note:: 泰勒公式(Taylor formula)、泰勒级数(Taylor series)、泰勒展开(Taylor explanation)、泰勒定理(Taylor theory)是一回事。 梯度下降法 ============================================ 我们把对数似然函数按照泰勒公式进行展开，但是我们只展开到一阶导数，把更高阶导数的和看做一个常数量constant， 则对数似然函数通过泰勒展开可以简化成： .. math:: f(x)_{Taylor} = f(x_0) + f'(x_0)(x-x_0) +constant 现在我们把对数似然函数按照泰勒展开公式进行操作： .. math:: :label: eq_34_30 \ell(\beta^{(t+1)}) = \ell(\beta^t) + \ell'(\beta^t)(\beta^{(t+1)} - \beta^t) + constant 其中，我们把 :math:`\beta^t` 看做是参数的初始值(或者说上一轮迭代后的值)，其值是已知的。 把 :math:`\beta^{(t+1)}` 看做是对数似然函数的极值点(或者说当前轮迭代后的值)，其值是未知的，是函数的自变量。 最大化对数似然函数是希望对数似然函数取得极大值点，那么就要求每一次迭代 :math:`\ell(\beta)` 的值要变得更大，直到达到极大值点。 那么就要满足： .. math:: \ell(\beta^{(t+1)}) - \ell(\beta^t) =\ell'(\beta^t)(\beta^{(t+1)} - \beta^t) +constant \ge 0 当等号成立的时候 :math:`\beta^{(t+1)}` 就是我们要找的极大值点，对上述公式进行移项处理，可得： .. math:: :label: eq_34_31 \beta^{(t+1)} = \beta^t - \frac{constant}{\ell'(\beta^t)} 通常我们把一阶导 :math:`\ell'(\beta^t)` 称为梯度(gradient)， :eq:`eq_34_31` 说明只要 :math:`\beta^{(t+1)}` 沿着 :math:`\beta^t` 的负梯度方向进行移动，我们终将能达到极值点。 注意 :math:`\frac{constant}{\ell'(\beta^t)}` **的绝对值的大小影响着前进的速度，** **其方向(正负号)决定目标函数是否向着极大值点移动。** 所以和下面的公式是等价的，:math:`\alpha` 称为学习率(learning rate)，是一个人工设置参数，控制的迭代的速度。 .. math:: :label: eq_34_32 \beta^{(t+1)} = \beta^t - \alpha \ell'(\beta^t) 利用 :eq:`eq_34_32` 进行参数迭代求解的方法就称为梯度下降法，梯度下降法的核心就是让参数(自变量)沿着负梯度的方向前进。 虽然理论上最终一定能到达极值点，但是实际上会受到学习率参数 :math:`\alpha` 的影响， 学习率可以理解成每次迭代前进的步长(step size)，步长越大前进的越快，收敛性速度就越快；反之，步长越小，收敛越慢。 但是步长如果大了，就会造成震荡现象，即一步迭代就越过了终点(极值点)，并且在极值点附近往返震荡，永远无法收敛。 为了保证算法能一定收敛，通常会为 :math:`\alpha` 设定一个较小的值。 关于 :math:`\alpha` 的更多讨论请参考其它资料。 .. _fg_34_3: .. figure:: pictures/34_3.png :scale: 70 % :align: center 梯度下降法中学习率的影响(图片来自网络) 牛顿法 ============================================ 梯度下降法虽然也能收敛到最优解，但是如果学习率设置(通常人工设置)不合理，可能会造成收敛速度太慢或者无法收敛的问题。 现在我们讨论另一中迭代算法，牛顿–拉夫森方法(Newton–Raphson)，一般简称牛顿法。 还是从泰勒展开公式开始，让我们考虑二阶泰勒展开： .. math:: :label: eq_34_33 \ell(\beta^{(t+1)}) = \ell(\beta^t) + \ell'(\beta^t)(\beta^{(t+1)} - \beta^t) + \frac{1}{2}\ell''(\beta^t)(\beta^{(t+1)} - \beta^t)^2 + constant 我们知道目标函数在极值点处的导数应该为0， 所以如果 :math:`\beta^{(t+1)}` 是极值点，那么有 :math:`\ell'(\beta^{(t+1)})=0` 。我们对 :eq:`eq_34_33` 进行求导，注意 :math:`theta_t` 和 :math:`\ell(\beta^t)` 都是已知常数量。 .. math:: \ell'(\beta^{(t+1)})= \ell'(\beta^t) + \ell''(\beta^t)(\beta^{(t+1)}-\beta^t)=0 通过移项可得： .. math:: :label: eq_34_34 \beta^{(t+1)} = \beta^t - \frac{\ell'(\beta^t)}{\ell''(\beta^t)} :math:`\ell''(\beta^t)` 是目标函数的二阶偏导数，由于参数 :math:`\beta` 是一个向量，所以二阶导是一个矩阵， 通常被称为海森矩阵(Hessian matrix)，常用符号 :math:`H` 表示。 .. math:: :label: eq_34_35 \beta^{(t+1)} = \beta^t - H(\beta^{(t)})^{-1} \ell'(\beta^t) 和梯度下降法的 :eq:`eq_34_32` 对比下发现，两者非常相似，不同的是牛顿法用Hessian矩阵的逆矩阵 :math:`H(\beta^{(t)})^{-1}` 替代了学习率参数，避免了需要人工设置学习率的问题。相比梯度下降法，牛顿法收敛速度更快，并且也没有震荡无法收敛的问题。 :math:`H(\beta^{(t)})` 是对数似然函数的二阶导数，可以在一阶导 :eq:`eq_34_jac` 的基础上再进行一次偏导得到。 .. math:: :label: eq_34_36 \left (\frac{\partial^2 \ell }{\partial \beta_j \partial \beta_k} \right ) &= \sum_{n=1}^N \frac{1}{a(\phi)} \left ( \frac{\partial}{\partial \beta_k} \right ) \left \{ \frac{y_n-\mu_n}{\nu(\mu_n)} \left ( \frac{\partial \mu}{\partial \eta} \right)_n x_{jn} \right \} &= \sum_{n=1}^N \frac{1}{a(\phi)} \left [ \left ( \frac{\partial \mu }{\partial \eta} \right )_n \left \{ \left ( \frac{\partial }{\partial \mu} \right )_n \left ( \frac{\partial \mu }{\partial \eta} \right )_n \left ( \frac{\partial \eta }{\partial \beta_k} \right )_n \right \} \frac{y_n-\mu_n}{\nu(\mu_n)} + \frac{y_n-\mu_n}{\nu(\mu_n)} \left \{ \left ( \frac{\partial }{\partial \eta} \right )_n \left ( \frac{\partial \eta }{\partial \beta_k} \right )_n \right \} \left ( \frac{\partial \mu }{\partial \eta} \right )_n \right ] x_{jn} &= -\sum_{n=1}^N \frac{1}{a(\phi)} \left [ \frac{1}{\nu(\mu_n)} \left ( \frac{\partial \mu}{\partial \eta} \right )_n^2 -(\mu_n-y_n) \left \{ \frac{1}{\nu(\mu_n)^2} \left ( \frac{\partial \mu }{\partial \eta} \right )_n^2 \frac{\partial \nu(\mu_n)}{\partial \mu} - \frac{1}{\nu(\mu_n)} \left ( \frac{\partial^2 \mu}{\partial \eta^2} \right )_n \right \} \right ] x_{jn}x_{kn} 有了一阶导数和二阶导数(Hessian matrix)，就可以按照 :eq:`eq_34_35` 迭代的更新参数值，直到算法收敛，最终得到参数的最大似然估计值。 完成优化后，必须为参数 :math:`\beta` 估计值估计一个合适的方差矩阵。 一个明显而合适的选择是基于估计的观测到的海森矩阵。 这是软件实现中最常见的默认选择，因为观察到的Hessian矩阵是估计算法的组成部分。 这样，它已经被计算并可用。 我们所介绍的 Newton–Raphson 算法并没有包括对 scale parameter， :math:`\phi` ，的估计功能，有兴趣的读者可以参考其它资料。 最终，Newton–Raphson 提供了如下功能： 1. 为所有单参数指数族的GLM成员模型提供一个参数估计算法。 2. 参数估计值的标准误(standard errors)的计算：负的 Hessian matrix 逆矩阵的对角线元素的平方根。 我们这里描述的 Newton–Raphson 算法不支持分散参数(dispersion parameter)， :math:`\phi` ，的估计，一些 Newton–Raphson 的扩展算法可以提供分散参数的ML估计。 然而，Newton–Raphson 方法存在两大缺点： 1. Hessian matrix 的计算复杂度非常高，:eq:`eq_34_36` 。 2. Hessian matrix 并不是一定存在逆矩阵的，当不存在逆矩阵时，算法无法执行。 由于这两个缺陷的存在，在实际应用中是很少直接使用 Newton–Raphson 法的， 而是使用牛顿法的改进算法。有很多牛顿法的改进算法被研究者提了出来， 多数改进算法的思路是找到一个海森矩阵的近似矩阵去替换原有的海森矩阵，进而避免上述两个缺点。 **迭代初始值的设定** 要实现 Newton–Raphson 迭代法， 我们必须对参数初始值有一个猜测。 但目前没有用于获得良好参数初值的全局机制， 但是有一个合理的解决方案可以在模型中存在"常数项系数"时获得起始值。 这里的"常数项"指的是线性预测器中截距部分 .. math:: \eta = \beta_0 \times 1 + \beta_1 x_1 +\dots + \beta_px_p 其中 :math:`\beta_0` 就是常数项系数。 如果模型包含常数项，则通常的做法是找到仅包含常数项系数的模型的估计值。 我们令： .. math:: \eta = \beta_0 然后令对数似然函数的一阶导数 :eq:`eq_34_jac` 为0，找到 :math:`\beta_0` 的解析解。 .. math:: :label: eq_34_37 \sum_{n=1}^N \frac{y_n-\mu_n}{a(\phi) \nu(\mu_n) } \left ( \frac{\partial \mu}{\partial \eta} \right )_n =0 通过上式是可以得到 :math:`\beta_0` 的一个估计值的。 比如如果是逻辑回归模型，则有： .. math:: a(\phi) &= 1 \nu(\mu) &= \mu(1-\mu) \mu &= \text{sigmoid}(\eta_n) = \text{sigmoid}(\beta_0) \frac{\partial \mu}{\partial \eta} &= \frac{\partial }{\partial \eta} \text{sigmoid} (\eta) = \mu(1-\mu) 代入到 :eq:`eq_34_37` 可得： .. math:: \sum_{n=1}^N \frac{(y_n- \mu_n ) }{\mu_n(1-\mu_n)} \mu_n(1-\mu_n) &= 0 &\Downarrow \sum_{n=1}^N (y_n- \mu_n) &=0 &\Downarrow \sum_{n=1}^N (y_n- \frac{1}{1+e^{-\beta_0}}) &=0 &\Downarrow \underbrace{\frac{1}{N}\sum_{n=1}^N y_n}_{\text{均值}\bar{y}} &= \frac{1}{1+e^{-\beta_0}} &\Downarrow{\text{sigmoid反函数求解}} \hat{\beta}_0 &= \ln \left ( \frac{\bar{y}}{1-\bar{y}} \right ) 然后我们就用 :math:`\beta=(\hat{\beta}_0,0,0,\dots,0)^T` 作为 Newton–Raphson 算法首次迭代时参数向量的初始值。 使用这种方法为我们提供了两个优点。 首先，我们从参数空间的合理子集开始迭代。 其次，对于ML，因为我们知道仅常数项系数模型的解决方案， 所以可以将训练模型的对数似然与算法初始步骤中获得的仅常数项系数模型的对数似然进行比较。 此比较是每个协变量（常数项系数除外）均为零的似然比检验，在下一节会介绍什么是似然比检验。 如果模型中没有常量项系数，或者如果我们无法解析法求解纯常数项系数模型，则必须使用更复杂的方法， 比如使用搜索方法寻找合理的初始点来开始我们的 Newton-Raphson 算法。 迭代重加权最小二乘(IRLS) ============================================ 我们已经知道在牛顿法中，参数是按照如下公式进行迭代更新的： .. math:: \beta^{(t+1)} = \beta^t - H(\beta^{(t)})^{-1} \ell'(\beta^t) 其中 :math:`H(\beta^{(t)})` 是对数似然函数的二阶偏导数， GLM的对数似然函数的二阶偏导数如 :eq:`eq_34_36` 所示， 显然其计算成本是高昂的，即使计算出来还面临着求逆矩阵的难题。 那么是不是可以找到 :math:`H(\beta^{(t)})^{-1}` 的一个近似替代值呢。 我们在 :numref:`ch_2_Fisher_Information` 已经讲过 Fisher information 的定义，我们知道信息矩阵(information matrix)就等于海森矩阵(Hessian matrix)的期望的负数。 .. math:: I(\beta) = - \mathop{\mathbb{E}}_{p(\mathcal{D} ; \beta)}[H(\beta^{(t)})] 我们可以用 :math:`I(\beta)` 代替 :math:`-H(\beta^{(t)})` ，回顾一下 :numref:`ch_2_Fisher_Information` 的内容， 信息矩阵 :math:`I(\beta)` 就是 Score function (对数似然函数的一阶导数) 的方差， 当然我们也可以直接对 :eq:`eq_34_36` 的求期望得到，结果是一样的。 .. math:: :label: eq_34_44 I(\beta) = - \mathbb{E}[H(\beta^{(t)})] &= \mathbb{E}[S(\beta)S(\beta)^T] &= \mathbb{E}[ \frac{ \partial \ell}{\partial \beta_j} \cdot \frac{ \partial \ell}{\partial \beta_k} ] &= \mathbb{E} \left [ \sum_{n=1}^N \frac{y_n-\mu_n}{a(\phi) \nu(\mu_n) } \left ( \frac{\partial \mu}{\partial \eta} \right )_n x_{jn} \cdot \sum_{n=1}^N \frac{y_n-\mu_n}{a(\phi) \nu(\mu_n) } \left ( \frac{\partial \mu}{\partial \eta} \right )_n x_{kn} \right ] &= \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } x_{jn} x_{kn} 现在我们把 :eq:`eq_34_35` 稍微变换一下，并且我们用 :math:`I(\beta^{(t)})` 替换 :math:`- H(\beta^{(t)})` 。 .. math:: :label: eq_34_45 \Delta \beta^{(t)} = \beta^{(t+1)} - \beta^{(t)} &= - H(\beta^{(t)})^{-1} \ell'(\beta^{(t)}) & \Downarrow{\text{移项}} - H(\beta^{(t)}) \Delta \beta^{(t)} &= \ell'(\beta^{(t)}) & \Downarrow{\text{Hessian的期望}} - \mathbb{E}[H(\beta^{(t)})] \Delta \beta^{(t)} &= \ell'(\beta^{(t)}) & \Downarrow{\text{信息矩阵}} I(\beta^{(t)})\Delta \beta^{(t)} &= \ell'(\beta^{(t)}) 现在我们把 :eq:`eq_34_jac` 和 :eq:`eq_34_44` 代入上式，得到一个等式： .. math:: :label: eq_34_46 \left \{ \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } x_{jn} x_{kn} \right \} \Delta \beta^{(t)} = \sum_{n=1}^N \frac{y_n-\mu_n}{a(\phi) \nu(\mu_n) } \left ( \frac{\partial \mu}{\partial \eta} \right )_n x_{n}^T \ \ \ \ \text{(等式A)} :eq:`eq_34_46` 先保留，我们记为等式A，我还需要借助另外一个等式。 我们知道，对于每条样本线性预测器的方程为： .. math:: \eta_n^{(t)} = x_n^T \beta^{(t)} 我们把 :eq:`eq_34_46` 中的 :math:`\Delta \beta^{(t)}` 换成 :math:`\beta^{(t)}` ，并结合线性预测器，可以得到如下等式。 .. math:: :label: eq_34_47 & \left \{ \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } x_{jn} x_{kn} \right \} \beta^{(t)} &= \left \{ \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } x_n x_n^T \right \} \beta^{(t)} &= \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } \beta^{(t)T} x_n x_n^T &= \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } \eta_n^{(t)} x_n^T 去掉 :eq:`eq_34_47` 的中间推导过程，直接得到如下等式，我们记为等式B。 .. math:: :label: eq_34_48 \left \{ \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } x_{jn} x_{kn} \right \} \beta^{(t)} = \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } \eta_n^{(t)} x_n^T \ \ \ \ \text{(等式B)} 现在我们把等式A( :eq:`eq_34_46` )和等式B( :eq:`eq_34_48` )的等号两边相加： .. math:: :label: eq_34_49 \left \{ \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } x_{jn} x_{kn} \right \} (\beta^{(t)} + \Delta \beta^{(t)} ) &= \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } \left \{ (y_n-\mu_n)\left ( \frac{\partial \eta}{\partial \mu}_n \right) + \eta_n^{(t)} \right \} x_n^T & \Downarrow \underbrace{ \left \{ \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } x_{jn} x_{kn} \right \} }_{\text{矩阵}} \beta^{(t+1)} &= \sum_{n=1}^N \left ( \frac{\partial \mu}{\partial \eta} \right )^2_n \frac{ 1}{ a(\phi) \nu(\mu_n) } \left \{ (y_n-\mu_n)\left ( \frac{\partial \eta}{\partial \mu} \right)_n+ \eta_n^{(t)} \right \} x_n^T 上式看上去很复杂，但其实可以转化成矩阵的乘法，我们定义如下两个矩阵： .. math:: W^{(t)} &= \text{diag} \left \{ \frac{ 1}{ a(\phi) \nu(\mu) } \left ( \frac{\partial \mu}{\partial \eta} \right )^2 \right \}_{(n\times n)} \ \ \ \ \text{对角矩阵} \mathcal{\eta}^{(t)} &= \left \{ (y-\mu) \left ( \frac{\partial \eta}{\partial \mu} \right) + \eta^{(t)} \right \}_{(n\times 1 )} 其中 :math:`\frac{\partial \mu}{\partial \eta}` 是激活函数的对 :math:`\eta` 的导数， :math:`\frac{\partial \eta}{\partial \mu}` 是链接函数对 :math:`\mu` 的导数。 :eq:`eq_34_49` 就可以改写成矩阵乘积的形式： .. math:: :label: eq_34_49_1 (X^TW^{(t)} X) \beta^{(t+1)} &= X^T W^{(t)} \eta^{(t)} & \Downarrow{\text{移项}} \beta^{(t+1)} &= (X^TW^{(t)} X)^{-1} X^T W^{(t)} \eta^{(t)} :eq:`eq_34_49_1` 就是最终参数向量更新的公式，它在形式上等价于加权的最小二乘法， 其中 :math:`W` 是权重矩阵，并且每一次迭代都要重新计算 :math:`W` ，所以我们把这个算法称为迭代重加权最小二乘法(Iteratively Reweighted Least Square,IRLS)， "reweighted" 指的就是每次迭代重新计算权重矩阵。 IRLS和牛顿法的区别就是，使用的 Hessian 矩阵的期望矩阵(expected Hessian)代替了原来的 Hessian 矩阵。 原来的 Hessian 矩阵是在观测样的基础上计算的，所以称为 observed Hessian ，对其在观测样本变量上求期望就得到了 expected Hessian 。 **迭代初始值的设定** 对比下 Newton–Raphson 算法的参数迭代公式( :eq:`eq_34_35` ) 和IRLS算法的参数迭代公式( :eq:`eq_34_49_1` )， 可以发现IRLS算法并不需要直接在 :math:`\beta^{(t)}` 的基础上进行参数迭代， IRLS算法的参数迭代仅仅依赖 :math:`\mu` 和 :math:`\eta` ，因此与 Newton–Raphson 算法不同的是，IRLS 不需要对参数向量 :math:`\beta` 进行初始值的猜测，只需要给 :math:`\mu` 和 :math:`\eta` 赋予一个初始值即可。 - 对于二项式分布，可以令 :math:`\mu_n^{(0)}=k_n(y_n+0.5)/(k_n+1)` ，:math:`\eta_n^{(0)}=g(\mu_n^{(0)})` 。 - 对于非二项式分布，可以令 :math:`\mu_n^{(0)}=y_n` ， :math:`\eta_n^{(0)}=g(\mu_n^{(0)})` 。 **收敛性判断** 在迭代的过程中，我们可以检查参数 :math:`\beta` 的相对变化来决定是否结束算法。 .. math:: \sqrt{\frac{ (\beta^{new}-\beta^{old})^T (\beta^{new}-\beta^{old}) }{ \beta^{old^T} \beta^{new} } } < \epsilon 也可以通过相对偏差来判断。 .. math:: \left|\frac{D(y-\mu^{new})-D(y,\mu^{old}) }{D(y,\mu^{old})} \right| <\epsilon **估计量的标准误差** 回顾一下在 :numref:`ch_2_MLE_estimator` 我们讲的最大似然估计量的评价， 我们知道最大似然估计量的协方差矩阵就是 :math:`I(\beta)^{-1}` ，显然在IRLS算法过程中已经计算出了 :math:`I(\beta)^{-1}=- \mathbb{E}[H(\beta)]=(X^TW^{(t)} X)^{-1}` ，所以使用IRLS我们可以很方便的得到估计量的标准误差。 .. math:: SSE = \sqrt{ \text{diag} [{(X^TW^{(t)} X)}^{-1} ]} **分散参数的估计** 虽然IRLS算法本身并不包含对分散参数的估计，但我们可以通过 Pearson :math:`\mathcal{X}^2` 统计来得到 :math:`a(\phi)` 。 .. math:: \hat{a}(\phi) = \frac{1}{N-p} \sum_{n=1}^N \frac{ (y_n - \hat{\mu}_n)^2}{\nu(\hat{\mu}_n)} 或者使用偏差： .. math:: \hat{a}(\phi) = \frac{D(y,\hat{\mu})}{N-p} :math:`N` 是观测样本的数量，:math:`p` 是参数向量 :math:`\beta` 的长度， :math:`\hat{\mu}_n` 是第n条样本的模型预测值， :math:`\hat{a}(\phi)` 服从自由度为 :math:`n-p` 的 :math:`\mathcal{X}^2_{n-p}` 分布。 goodness of fit ################################################# 我们知道，模型的参数越多对数据的拟合程度就越好，极端情况下，模型参数的数量和样本的数量相同， 这时就相当于对每条样本都有一个独立的参数(模型)去拟合它，理论上可以完美拟合所有的样本。 我们把这样的模型成为之饱和模型(saturated model)，也可以称为 完整模型(full model)或者最大模型(maximal model)。 饱和模型虽然能完美拟合数据集，但它并没有从数据集中学习出任何的统计信息(统计规律)，所不具备泛化能力， 俗称过拟合(over-fitted)。 通过为饱和模型中的参数添加约束，比如令一些参数值为0，相当于去掉了一个参数，这样就得到了简化的模型。 简化模型对数据集拟合度下降了，但是其泛化能力会得到提升， 更少的参数数量可以得到更大的泛化能力。 但是参数数量变少，会降低拟合程度，参数数量越少拟合度就越差，所以也不是参数越少越好。 在开发一个模型时，我们希望模型的预测值 :math:`\hat{y}` 尽可能的接近数据的真实值 :math:`y` ，对于一个规模为N的观测值样本，我们可以考虑参数数量在 :math:`[1,N]` 之间的候选模型， 最简单的模型是只有一个参数的模型，但它对所有的样本的预测值都是一样的，缺乏拟合能力。 最复杂的模型是含有N个参数的模型，它可以完美拟合所有样本，但是它缺乏泛化能力。 理论上，我们期望得到一个参数数量尽可能少，又能保持拟合能力的模型。 在统计学中，似然比检验(likelihood-ratio test,LRT)用来对比两个模型对于当前数据集的拟合程度， 其是利用似然比(likelihood-ratio LR)的值来比较两个模型的优劣。 LRT的计算公式如下： .. math:: LR = 2 \ln \frac{L1}{L2} = 2 (\ln L1-\ln L2) 其中L1为复杂模型最大似然值，L2为简单模型最大似然值。 从公式可以看出，似然比就是两个模型的似然值之比的对数，也可以看成是两个模型对数似然值的差值。 似然(likelihood)，实际上也可以翻译为可能性，表示的是样本发生的概率，显然似然值越大的模型对数据的拟合也就越好。 似然比就是直接比较两个模型的似然值大小。 但是并不是任意两个模型都可以应用似然比去比较，只有在特定条件下似然比才有意义。 1. 两个模型采用同一份数据集，样本的数量和特征都是相同的。这很好理解，不同数据集似然值自然是不同的，没有比较的意义。 2. 两个模型是嵌套关系(nested)。所谓嵌套关系就是，其中一个模型是通过把另一个模型中的部分参数设置为0而得到的。 当样本足够大时，似然比是渐进服从卡方分布的，其自由度等于两个模型的参数数量之差(参数值为0的参数的数量)。 这样根据卡方分布临界值表，我们就可以判断模型差异是否显著。 在GLM中，我们定义一个评估模型拟合优度(Goodness of fit)的指标，称之为偏差(deviance)。 **偏差的计算方法就是饱和模型和拟合模型之间的似然比。** 我们用符号 :math:`L_m` 表示我们拟合出的目标模型的似然值， 用符号 :math:`L_f` 表示饱和模型的似然值。 用符号 :math:`D` 表示偏差，其通过下式给出： .. math:: :label: eq_34_50 D = 2 (\ln L_f -\ln L_m) 我们知道在GLM中，模型的预测值 :math:`\hat{y}` 就是分布 :math:`p(y|x)` 的期望值 :math:`\mathbb{E}[p(y|x)]=\hat{\mu}` ，即 :math:`\hat{y}=\hat{\mu}` 。 所以这里我们用 :math:`\hat{\mu}` 表示模型的预测值。 并且我们知道，指数族形式的规范参数 :math:`\theta` 可以看做是一个关于均值参数 :math:`\mu` 的函数， 因此，在GLM框架下，目标模型的似然值为： .. math:: :label: eq_34_51 L_m = \exp \left \{ \sum_{n=1}^N \frac{y_n \theta(\hat{\mu}_n) -b(\theta(\hat{\mu}_n))}{a(\phi)} + \sum_{n=1}^N c(y_n;\phi) \right \} 对于饱和模型，每条样本的预测值就是样本的真实值，即 :math:`\hat{y_n}=y_n` ，换句话说， 对于饱和模型，满足 :math:`\hat{y_n}=\hat{\mu}_n=y_n` 。因此，饱和模型的似然值为： .. math:: :label: eq_34_52 L_f = \exp \left \{ \sum_{n=1}^N \frac{y_n \theta(y_n) -b(\theta(y_n))}{a(\phi)} + \sum_{n=1}^N c(y_n;\phi) \right \} 把 :eq:`eq_34_51` 和 :eq:`eq_34_52` 代入到 :eq:`eq_34_50` 可得到GLM的偏差： .. math:: D &= 2 (\ln L_f -\ln L_m) &= \frac{2}{a(\phi)} \sum_{n=1}^N [ y_n \{ \theta(y_n) - \theta(\hat{\mu}_n) \} - b\{\theta(y_n)\} + b\{\theta(\hat{\mu}_n)\} ] &\triangleq 2 \sum_{n=1}^N [ y_n \{ \theta(y_n) - \theta(\hat{\mu}_n) \} - b\{\theta(y_n)\} + b\{\theta(\hat{\mu}_n)\} ] 由于 :math:`a(\phi)` 与样本和参数都无关，是一个常量，所以通常会被省略掉。 对于一个特定的样本集和模型，饱和模型的似然值是一个常数值， 因此 **最大化似然函数和最小化偏差是等价的** **，在进行参数学习时可以用最小化偏差替代最大化对数似然函数。** 在统计中，偏差(deviance)是统计模型的拟合优度统计； 它通常用于统计假设检验。 它是将普通最小二乘中的残差平方和用于通过最大似然实现模型拟合的情况的概括。 它在指数弥散模型和广义线性模型中起着重要作用。 https://newonlinecourses.science.psu.edu/stat504/node/220/ https://en.wikipedia.org/wiki/Deviance_(statistics) https://www.xiaofandajie.top/2018/03/06/%E6%9E%81%E5%A4%A7%E4%BC%BC%E7%84%B6%E6%B3%95%E4%B8%8Elogistic%E5%9B%9E%E5%BD%92/ https://bookdown.org/hezhijian/book/est.html#section-30 连续值响应模型 ################################################# 我们知道在大部分机器学习问题中，都是在寻找因变量 :math:`Y` 和自变量 :math:`X` 之间的关系， 在概率的语义下，用条件概率 :math:`P(Y|X)` 定义这种关系。 GLM框架给了定义条件概率 :math:`P(Y|X)` 的一种通用性方法， GLM包含了一类统计模型，可以在不同的场景下选择其中合适的模型去应用。 在GLM框架下要确定一个具体的模型，理论上需要确定两个信息： - 根据标签 :math:`Y` 的数据分布选取一个合适的指数族分布作为变量 :math:`Y` 的概率分布假设。 - 确定一个连接函数 :math:`g(\cdot)`，把特征数据 :math:`X` 的线性预测器 :math:`\beta^T x` 与 :math:`Y` 的概率分布的均值参数 :math:`\mu` 连接在一起 :math:`\beta^T x=g(\mu)` 。 不同的场景标签 :math:`Y` 拥有不同的数据范围和分布，就需要选取特定的指数分布。 本节开始，我们介绍指数族中常见分布的GLM，帮助大家在遇到具体的场景时， 能用GLM解决问题。根据数据的不同，我们分为如下几类： - 连续值变量，对应着实数值回归问题场景。 - 二值离散变量，对应着二分类问题场景。 - 多值离散变量，对应着多分类问题场景。 - 计数离散变量。 高斯族 ===================================================== 1800年代， 约翰·卡尔·弗里德里希·高斯（Johann Carl Friedrich Gauss） 描述了最小二乘拟合方法和以他的名字命名的分布。 高斯密度函数具有对称的钟形形状，通常称为正态分布。 正态累积分布函数是指数分布族的成员，因此可以用于GLM框架。 此外，由于GLM的理论最初被认为是对普通最小二乘（OLS）模型的扩展， 因此我们将首先说明OLS如何适合GLM框架。 这样一来，了解其他GLM模型如何概括此基本形式就变得容易了。 高斯分布是实数域的连续分布，其输入域是整个实数域 :math:`R=(-\infty,+\infty)` 。用均值 :math:`\mu` 进行参数化的高斯概率密度函数可以表示为： .. math:: f(y;\mu,\sigma^2)=\frac{1}{\sqrt{2\pi \sigma^2}} \exp \left \{ - \frac{(y-\mu)^2}{2\sigma^2} \right \} 其中 :math:`f(\cdot)` 表示在给定参数 :math:`\mu` 和 :math:`\sigma^2` 时，变量 :math:`y` 的概率密度函数的一般形式。 :math:`y` 表示输出变量， :math:`\mu` 表示均值参数， :math:`\sigma^2` 表示 scale parameter 。 基于高斯或正态分布的回归模型通常称为OLS模型，该标准回归模型通常是统计学模型的入门模型。 在高斯回归模型中，我们定义响应变量和特征变量的之间的关系是： .. math:: y = \beta^T x + \epsilon 其中 :math:`\beta^T x` 是输入变量的线性预测器，:math:`\epsilon \sim \mathcal{N}(0,\sigma^2)` 是一个高斯误差项， 输出变量 :math:`y` 是一个服从高斯分布的变量 :math:`y \sim \mathcal{N}(\beta^Tx,\sigma^2)` 。现在我们用GLM框架来解释OLS。 GLM中高斯族的定义 ----------------------------------------------------------------------- 当响应变量 :math:`Y` 数值范围是实数值时，并且其数据分布是近似高斯分布时， 我们就可以为响应变量 :math:`Y` 建立高斯分布假设。 在GLM框架下，我们需要用指数分散族(EDF)的形式表示指数族的分布，高斯分布的EDF表示为： .. math:: f(y;\mu,\sigma^2) &=\frac{1}{\sqrt{2\pi \sigma^2}} \exp \left \{ - \frac{(y-\mu)^2}{2\sigma^2} \right \} &= \exp \left \{ -\frac{(y-\mu)^2}{2\sigma^2} - \frac{1}{2} \ln (2\pi\sigma^2) \right \} &= \exp \left \{ -\frac{y\mu - \mu^2/2}{\sigma^2} - \frac{y^2}{2\sigma^2} - \frac{1}{2} \ln (2\pi\sigma^2) \right \} 和EDF的标准形式 :eq:`eq_34_EDF` 对比下，有： .. math:: \theta &= \mu b(\theta) &= \mu^2/2 a(\phi) &= \sigma^2 对于高斯分布的GLM来说，链接函数 :math:`g` 是恒等函数，即 :math:`\eta=\mu` ，我们称之为规范链接函数(canonical link)。 此外，对于高斯分布，规范参数 :math:`\theta` 和均值参数 :math:`\mu` 之间也是恒等函数关系， 所以，对于高斯分布有： .. math:: \beta^Tx = \eta=\mu=\theta 均值参数可以通过累积函数 :math:`b(\theta)` 的一阶导数求得： .. math:: b'(\theta) &= \frac{\partial b}{\partial \theta} &= \frac{\partial b}{\partial \mu} \frac{\partial \mu}{\partial \theta} &= (\mu)(1) = \mu 方差可以通过累积函数 :math:`b(\theta)` 的二阶导数和 :math:`a(\phi)` 的乘积求得： .. math:: b''(\theta) &= \frac{\partial^2 b}{\partial \theta^2} &= \frac{\partial }{\partial \theta} \left ( \frac{\partial b}{\partial \mu} \frac{\partial \mu}{\partial \theta} \right ) &= \frac{\partial }{\partial \theta}(\mu) &= \frac{\partial }{\partial \mu} \mu \frac{\partial \mu}{\partial \theta} &= (1)(1)=1 高斯分布的模型预测值为： .. math:: \hat{y} = \mathbb{E}[p(y|x)] = \mu =\beta^Tx 这和OLS是等价的。 现在我们来看下高斯分布的GLM模型的偏差的计算。 观测样本数据的对数似然函数可以通过下式计算： .. math:: \ell(\hat{\mu},\sigma^2;y)=\sum_{n=1}^N \left \{ \frac{y_n \hat{\mu}_n - \hat{\mu}_n^2/2}{\sigma^2} - \frac{y_n^2}{2\sigma^2} - \frac{1}{2} \ln (2\pi\sigma^2) \right \} 饱和模型的对数似然函数为； .. math:: \ell(y,\sigma^2;y) &=\sum_{n=1}^N \left \{ \frac{y_n^2 - y_n^2/2}{\sigma^2} - \frac{y_n^2}{2\sigma^2} - \frac{1}{2} \ln (2\pi\sigma^2) \right \} &= \sum_{n=1}^N \left \{- \frac{1}{2} \ln (2\pi\sigma^2) \right \} 模型的偏差为： .. math:: D &= 2\{ \ell(y,\sigma^2;y) - \ell(\hat{\mu},\sigma^2;y) \} &= 2 \sum_{n=1}^N \left \{ - \frac{1}{2} \ln (2\pi\sigma^2) - \frac{y_n \hat{\mu}_n - \hat{\mu}_n^2/2}{\sigma^2} + \frac{y_n^2}{2\sigma^2} + \frac{1}{2} \ln (2\pi\sigma^2) \right \} &= 2 \sum_{n=1}^N \left \{\frac{ y_n^2 - 2y_n \hat{\mu}_n + \hat{\mu}_n^2 }{2\sigma^2} \right \} &= \frac{2}{2\sigma^2} \sum_{n=1}^N \left \{ y_n^2 - 2y_n \hat{\mu}_n + \hat{\mu}_n^2 \right \} &\triangleq \sum_{n=1}^N (y_n-\hat{\mu}_n)^2 我们发现高斯偏差等价于平方损失的和，这是高斯族独有的属性。 ML 估计 -------------------------------------- 现在我们来看下高斯分布的GLM的似然估计，首先我们需要把均值参数化的模型转化成使用 :math:`\beta^Tx` 参数化的形式。 .. math:: f(y;\mu,\sigma^2) &=\frac{1}{\sqrt{2\pi \sigma^2}} \exp \left \{ - \frac{(y-\mu)^2}{2\sigma^2} \right \} &= \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left \{ - \frac{(y-\beta^Tx)^2}{2\sigma^2} \right \} 对数似然函数为： .. math:: \ell(\beta,\sigma^2;y) = \sum_{n=1}^N \left \{ \frac{y_n \beta^Tx_n -(\beta^Tx_n)^2/2}{\sigma^2} -\frac{y_n^2}{2\sigma^2} - \frac{1}{2} \ln(2\pi\sigma^2) \right \} 对于采用规范链接函数的高斯模型的对数似然函数的一阶偏导数为： .. math:: \frac{\partial \ell}{\partial \beta_j} &= \sum_{n=1}^N \frac{1}{\sigma^2}(y_n-\beta^T x_n)x_{jn} \frac{\partial \ell}{\partial \sigma} &= \sum_{n=1}^N \frac{1}{\sigma} \left \{ \left ( \frac{y_n-\beta^T x_n}{\sigma} \right )^2 -1 \right \} 二阶偏导数为： .. math:: \frac{\partial^2 \ell}{\partial \beta_j\beta_k} &= - \sum_{n=1}^N \frac{1}{\sigma^2}x_{jn}x_{kn} \frac{\partial^2 \ell}{\partial \beta_j \partial \sigma} &= -\sum_{n=1}^N \frac{2}{\sigma^3} (y_n-\beta^Tx_n) x_{jn} \frac{\partial^2 \ell}{\partial \sigma \partial \sigma} &= -\sum_{n=1}^N \frac{1}{\sigma^2} \left \{ 3\left ( \frac{y_n-\beta^T x_n}{\sigma} \right )^2 - 1\right \} 对数高斯(log-Gaussian)模型 -------------------------------------------- 使用GLM作为模型构建框架的重要原因是能够轻松调整模型以适合特定的标签数据情况。 规范链接(canonical-link)高斯模型假定标签为正态分布。 尽管正态分布模型对于克服此假设具有一定的鲁棒性，但仍然有许多数据情况不适合正态分布模型的情况。 不幸的是，许多研究人员已将规范链接高斯模型用于不符合高斯模型所基于假设的数据情况。 当然，许多研究人员很少接受非正态分布模型建模方面的培训。 现在，大多数流行的软件包都具有GLM功能，或者至少实现了许多最广泛使用的GLM程序， 包括 logistic, probit, and Poisson regression。 线性预测器 :math:`\eta=\beta^Tx` 的取值范围是整个实数域 :math:`R=(-\infty,+\infty)` ，链接函数的主要作用就是将实数域的 :math:`\eta` 映射到标签数据的值域。 当标签数据仍然服从高斯分布，但是其范围不再是整个实数域，而是大于0的实数域时，恒等函数的规范链接(canonical-link)函数将不再适用， 这时就需要选取一个合适的链接函数，将 :math:`\eta` 从 :math:`R=(-\infty,+\infty)` 映射到 :math:`R=(0,+\infty)` 。 log-Gaussian模型仍然是基于高斯分布，但它的链接函数不再是规范链接(canonical-link)， 而是对数链接函数，对数链接通常用于标签数据是非负值的情况。 .. math:: \eta &= \ln \mu \mu &= \exp \{ \eta\} 在提出GLM框架之前，研究人员在遇到标签数据都是正数的情形时，采用的方法是把标签数据进行对数转化， :math:`y=\ln(y)` ，通过这种转换令数据符合整个实数域上的高斯模型假设，然后再利用标准高斯模型进行建模。 然而事实证明，这种方法在易用性和效果上都不如采用对数链接函数的对数高斯模型。 对数高斯模型的实现非常简单，只需要把标准高斯模型（采用恒等函数，规范链接函数）的链接函数替换成对数链接函数即可。 .. math:: f(y;\beta,\sigma^2) &=\frac{1}{\sqrt{2\pi \sigma^2}} \exp \left \{ - \frac{(y-\exp(\beta^Tx))^2}{2\sigma^2} \right \} &= \exp \left \{ -\frac{(y-\exp(\beta^Tx))^2}{2\sigma^2} - \frac{1}{2} \ln (2\pi\sigma^2) \right \} &= \exp \left \{ -\frac{y\exp(\beta^Tx) - \exp(\beta^Tx)^2/2}{\sigma^2} - \frac{y^2}{2\sigma^2} - \frac{1}{2} \ln (2\pi\sigma^2) \right \} 对数高斯模型的对数似然函数(log-likelihood)为： .. math:: \ell(\mu;y,\sigma^2) = \sum_{n=1}^N \left [ \frac{y_n \exp(\beta^Tx_n) - \{\exp(\beta^Tx_n)\}^2/2 }{\sigma^2} -\frac{y_n^2}{2\sigma^2} - \frac{1}{2}\ln(2\pi\sigma^2) \right ] 对数高斯模型就是在标准高斯模型的基础上，把链接函数由 :math:`\eta=\mu` 改成了 :math:`\eta=\ln(\mu)` ，相应的反链接函数由 :math:`\mu=\eta` 变成 :math:`\mu=\exp(\eta)` 。 规范链接函数的导数为1，然而对数链接函数的导数为 :math:`1/\mu` 。 注意对 :eq:`eq_34_jac` 相应位置替换即可得到似然函数的梯度。 Gamma族 ===================================================== Gamma函数 --------------------------------------------------- Gamma模型用于对输出数据只能是大于等于0的情况， 二项响应模型 ################################################# 多项响应模型 ################################################# 计数响应模型 ################################################# 泊松分布就是描述某段时间内，事件具体的发生概率。 泊松分布 ============================================== 泊松分布实际上是二项分布的试验(trials)次数n趋近于无穷时的场景，我们用一个例子说明。 假设一个交通观察员需要对某个路口的车流量进行建模，然后用模型预测未来一个小时内从这个路口通过的车次。 为了简化问题，我们假设路口的交通量不存在高峰期低峰期，即交通量不会随着时间的变化而变化， 并且每个时间片段内通过的车辆是互不影响的，当前一小时内车辆通过与否不影响下一个小时车辆。 观察员首先根据这个路口历史上车辆通过情况，计算出平均每小时通过车辆的数量为 :math:`\lambda` 。我们把一个小时内从路口通过的车次数看做一个随机变量，用符号 :math:`X` 表示， 那么 :math:`\lambda` 就是变量 :math:`X` 的数学期望。 .. math:: \mathbb{E}[X] = \lambda 我们把一辆车通过与否看做是一个伯努利变量，类似投硬币实验，1表示车辆通过，0表示车辆不通过。 把一个小时的时间区间均分成n个时间片段，比如每分钟作为一个片段，这时 :math:`n=60` 。 每个时间片段有车辆通过就是一次成功的实验(类似于投硬币正面向上)， 没有车辆通过就是一次失败的实验(类似于投硬币反面向上)， 这样就把一小时内车辆通过问题转化成一个二项分布问题， 在n次实验中有k次成功(车辆通过)概率分布函数可以写成： .. math:: :label: eq_34_201 p(X=k) = \binom{n}{k} p^k \left ( 1-p \right )^{n-k} 其中 :math:`p` 是一次实验的成功概率，我们已经通过历史数据知道平均每小时(n次实验)中通过的车次数为 :math:`\lambda` ，意味着n次实验中有 :math:`\lambda` 次成功， 单次实验成功的概率(平均一分钟内通过车辆数)为： .. math:: p = \frac{\lambda}{n} 但是我们并不能保证每分钟只有一辆车通过，我们需要保证一个时间片段内只有一辆车通过(一次实验) 以上的二项分布的假设才有意义。 理论上，我们只要把一小时的时间区间拆的足够小，比如拆成每秒，甚至是每毫秒为一个时间片段， 这样就能尽量保证每个时间片段内只会有一辆车通过。 :math:`n` 越大时间片段就越小，极限情况，我们可以把一小时分割成每个车辆通过的"瞬间"。 换句话说，只要 :math:`n \to \infty` 上述假设就是成立的，因此我们为 :eq:`eq_34_201` 加上极限操作。 .. math:: :label: eq_34_202 p(X=k) = \lim_{n \to \infty} \binom{n}{k} p^k \left ( 1-p \right )^{n-k} 我们发现 :eq:`eq_34_202` 就是二项分布的极限情况，表示的是路口未来一小时内通过的车辆数的概率分布， :math:`p(X=k)` 表示在一小时内通过车辆数为 :math:`k` 的概率。 :math:`\lambda` 表示这个时间区间内通过车辆数的期望值， 至于这个时间区间是一小时还是两小时并不重要， 只要是一个固定的时间区间就行， 所以可以看成是单位时间区间内，或者 :math:`t` 时间区间内。 :eq:`eq_34_202` 带有极限操作，事实上可以通过一些变换去掉极限符号， 现在我们尝试对其进行一些变换。 .. math:: :label: eq_34_210 p(X=k) &= \lim_{n \to \infty} \binom{n}{k} p^k \left ( 1-p \right )^{n-k} &= \lim_{n \to \infty} \frac{n!}{(n-k)!k!} \left (\frac{\lambda}{n} \right )^k \left ( 1-\frac{\lambda}{n} \right )^{n-k} &= \lim_{n \to \infty} \frac{n!}{(n-k)!k!} \frac{\lambda^k}{n^k} \left ( 1-\frac{\lambda}{n} \right )^n \left ( 1-\frac{\lambda}{n} \right )^{-k} .. math:: :label: eq_34_211 \frac{n!}{(n-k)!} = \frac{n (n-1) \cdots 2\times 1}{(n-k)(n-k-1) \cdots 2\times 1} = \underbrace{n (n-1) \cdots (n-k+1)}_{\text{k个}} .. math:: :label: eq_34_212 \lim_{x \to a} f(x)g(x) = \lim_{x \to a} f(x) \lim_{x \to a}g(x) :eq:`eq_34_210` 变成： .. math:: :label: eq_34_213 p(X=k) &= \lim_{n \to \infty} \frac{n (n-1) \cdots (n-k+1)}{n^k} \frac{\lambda^k}{k!} \left ( 1-\frac{\lambda}{n} \right )^n \left ( 1-\frac{\lambda}{n} \right )^{-k} &=\frac{\lambda^k}{k!} \lim_{n \to \infty} \left [ \frac{n (n-1) \cdots (n-k+1)}{n^k} \right ] \lim_{n \to \infty} \left [ \left ( 1-\frac{\lambda}{n} \right )^n \right ] \lim_{n \to \infty}\left [ \left ( 1-\frac{\lambda}{n} \right )^{-k}\right ] .. math:: :label: eq_34_214 \lim_{n \to \infty} \frac{n (n-1) \cdots (n-k+1)}{n^k} = 1 \lim_{n \to \infty} \left ( 1-\frac{\lambda}{n} \right )^n = e^{\lambda} \lim_{n \to \infty}\left [ \left ( 1-\frac{\lambda}{n} \right )^{-k}\right ] = 1 .. math:: :label: eq_34_215 p(X=k|n) = \frac{\lambda^k}{k!} \times 1 \times e^{\lambda} \times 1 =\frac{\lambda^k}{k!} e^{\lambda} 各个时间区间之间是相互独立的，互不影响，也就是不会因为当前时间区间内有车辆通过，而导致下一个时间区间内通过的车辆受到影响。 泊松分布的应用并不是仅限于固定的时间区间，理论上只要是固定的区间(fixed interval)即可， 比如固定大小的时间、长度、空间、面积、体积等等。 GLM扩展 #################################################