广义线性模型（GLMs）及算法介绍

一般我们了解的线性模型是针对连续性变量，并且服从正态分布的，但是在实际应用上显得非常的局限。因为我们我看到的数据很多都是离散的，而且不是服从正态分布的。针对这种情况，对传统线性模型进行推广，行成了现在的广义线性模型。广义线性模型使得变量从正态分布拓展到指数分布族，从连续型变量拓展到离散型变量，这就使得在现实中有着很好的运用，下面开始介绍广义线性模型。###广义线性模型（GLM）定义由以下三..

阴天了

17281人浏览 · 2019-02-27 17:13:04

阴天了 · 2019-02-27 17:13:04 发布

广义线性模型（GLM）定义

由以下三部分组成：
1 随机部分
随机样本 $Y_{1},Y_{2},...,Y_{n}$ 服从的分布来自指数分布族，即 $Y_{i}$ 的分布形式为 $f(yi;θi,ϕ)=exp{yiθi−b(θi)a(ϕ)+c(yi,ϕ)}f\left ( y_{i};\theta _{i},\phi \right )=exp\left \{ \frac{y_{i}\theta _{i}-b\left ( \theta _{i} \right )}{a\left ( \phi \right )}+c\left ( y_{i},\phi \right ) \right \}$
其中参数 $θi\theta_{i}$ 叫做正则参数，并且随着指数 $i （ i = 1, 2, . . ., n ）$ 而变化，但是扰乱因子 $ϕ\phi$ 是个常数。
2 系统部分
对于第 $i$ 个观测 $Y_{i}$ ，我们有一个称为系统部分的线性预测值，即所研究变量的线性组合，即 $ηi=xiTβ=∑j=1pxijβj,i=1,2,...n\eta _{i}=x_{i}^{T}\beta =\sum_{j=1}^{p}x_{ij}\beta _{j},i=1,2,...n$
3 连接函数
有一个单调可微函数 $g()g\left ( \right )$ 称为连接函数，它将随机部分的期望和系统部分连接起来，通过下面的等式 $g(μi)=ηi=xiTβ,i=1,2,...n,g\left ( \mu_{i} \right )=\eta _{i}=x_{i}^{T}\beta ,i=1,2,...n,$ 其中 $μi=E(Yi)\mu_{i}=E\left ( Y_{i} \right )$ 是 $Y_{i}$ 的期望。
矩阵表示：
$η=[η1η2⋮ηn]n×1,μ=[μ1μ2⋮μn]n×1,X=[x1′x2′⋮xn′]n×p\eta =\begin{bmatrix}\eta _{1}\\ \eta _{2}\\ \vdots \\ \eta _{n}\end{bmatrix}_{n\times 1},\mu=\begin{bmatrix}\mu _{1}\\ \mu _{2}\\ \vdots \\ \mu _{n}\end{bmatrix}_{n\times 1},X=\begin{bmatrix}x_{1}^{'}\\ x_{2}^{'}\\ \vdots \\ x_{n}^{'}\end{bmatrix}_{n\times p}$
那么连接函数可以用矩阵形式表示 $g(μ)=η=Xβg\left ( \mu \right )=\eta =X\beta$

连接函数介绍

1、正如 $G L M s$ 的定义所指出的那样，连接函数是单调可微的，用于连接随机部分的期望和系统部分的线性预测值
2、选择与分布相关的自然参数作为连接函数，在这种情况下，它被称为点则连接。具体而言，如果连接函数 $g ()$ 采用与自然参数相同的函数形式，那么它被称为点则连接函数。
3、点则连接的优点是它可以带来非常好的统计特性，并且使用起来很方便。例如，对于最常用的分布，我们有以下点则连接函数。

Normal	$μ=η\mu =\eta$	(identity-link)
Poisson	$logμ=ηlog\mu =\eta$	(log-link)
Bernoulli	$logπ1−π=ηlog\frac{\pi}{1-\pi}=\eta$	(logistic-link)
Binomial	$logπ1−π=ηlog\frac{\pi}{1-\pi}=\eta$	(logistic-link)

4、然而，点则连接并不是连接函数的唯一选择。其他可能的 $G L M s$ 连接函数包括
（a）二项分布的Probit连接： $η=Φ−1(π)\eta =\Phi ^{-1}\left ( \pi \right )$ ; $0<π<10<\pi<1$ ,其中 $Φ()\Phi()$ 叫做累计分布函数（不是点则连接呦）
（b）补充的二项分布的log-log连接 $η=log{−log(1−π)},0<π<1\eta =log\left \{ -log\left ( 1-\pi \right ) \right \},0<\pi<1$
（c）属于任何幂族分布的连接 $η={μλ,ifλ≠0logμ,ifλ=0\eta =\left\{\begin{matrix}\mu ^{\lambda }, if \lambda \neq 0 & \\ log\mu ,if \lambda =0&\end{matrix}\right.$

最大似然估计（MLE）的一般原则

假设我们对未知参数 $θ\theta$ 的对数似然函数，比如说 $l(θ)l\left ( \theta \right )$ 我们想找出 $θ\theta$ 的最大似然估计（MLE） $θ^\hat{\theta }$ ，即 $θ^≡argmaxθ⊂Ω{l(θ)}\hat{\theta }\equiv arg \underset{\theta \subset \Omega }{max}\left \{ l\left ( \theta \right ) \right \}$
这是估计方程的解 $∂l(θ)∂θ=0\frac{\partial l\left ( \theta \right )}{\partial \theta }=0$
1、在这种情况下，例如，对于正态分布参数 $θ\theta$ 的最大似然估计 $θ^\hat{\theta }$ 可以有一个明确的数学表达式 $（μ=1n∑i=1nlnxi）（\mu =\frac{1}{n}\sum_{i=1}^{n}lnx_{i}）$
2、一般来说，最大似然估计 $θ\theta$ 没有没有明确的数学解。相反，需要某些数值优化方法。
3、统计学中最常用的两种优化方法是Newton-Raphson算法和Fisher得分算法，他们都涉及计算 $l(θ)l\left ( \theta \right )$ 对 $θ\theta$ 的2次偏导数。

Newton-Raphson算法

该算法计算最大似然估计 $θ^\hat{\theta }$ ，通过下面的迭代：
$θm=θm−1+[−l′′(θ(m−1))]−1[l′(θ(m−1))]（1）\theta ^{m}=\theta ^{m-1}+\left [ -l^{''} \left ( \theta ^{\left ( m-1 \right )} \right )\right ]^{-1}\left [ l^{'} \left ( \theta ^{\left ( m-1 \right )} \right )\right ]（1）$
其中 $m = 1, 2, . . .$ 这里， $l′(θ(m−1))=∂l(θ)∂θ∣θ=θ(m−1)l^{'}\left ( \theta ^{\left ( m-1 \right )} \right )=\frac{\partial l\left ( \theta \right )}{\partial \theta }|_{\theta =\theta ^{\left ( m-1 \right )}}$ $l′′(θ(m−1))=∂2l(θ)∂θ∂θT∣θ=θ(m−1)l^{''}\left ( \theta ^{\left ( m-1 \right )} \right )=\frac{\partial ^{2}l\left ( \theta \right )}{\partial \theta \partial \theta ^{T}}|_{\theta =\theta ^{\left ( m-1 \right )}}$ 是 $p×1p\times 1$ 和 $p×pp\times p$ 的向量和矩阵 $（p是θ的维数）（p是\theta的维数）$
注1 $l′(θ)被称为θ的得分函数。−l′′(θ)被称为θ的观测信息矩阵l^{'}\left ( \theta \right )被称为\theta的得分函数。-l^{''}\left ( \theta \right )被称为\theta的观测信息矩阵$
注2 算法（1）需要初始值，例如 $θ0\theta ^{0}$ ，以开始迭代过程。初始值的选择通常需要经验。
注3 算法（1）迭代直到收敛。例如，当迭代结果满足 $∥θ(m)−θ(m−1)∥∥θ(m−1)∥≤10−5\frac{\left \| \theta ^{\left ( m \right )}-\theta ^{\left ( m-1 \right )} \right \|}{\left \| \theta ^{\left ( m-1 \right )}\right \|}\leq 10^{-5}$ 则迭代停止。 $θ(m)\theta ^{\left ( m \right )}$ 可以认为是最大似然估计 $θ^\hat{\theta }$ 。

Fisher得分算法

Fisher得分算法与Newton-Raphson算法相同，只是（1）式中的观测矩阵 $−l′′(θ)-l^{''}\left ( \theta \right )$ 被Fisher信息矩阵所代替 $I(θ)=E[−l′′(θ)]=−∫l′′(θ∣Y)fY(Y∣θ)dYI\left ( \theta \right )=E\left [ -l^{''}\left ( \theta \right ) \right ]=-\int l^{''}\left ( \theta |Y \right )f_{Y}\left ( Y|\theta \right )dY$
注释 Fisher得分算法和Newton-Raphson算法一般收敛于同一解。对于前者，在某些情况下，在信息矩阵的解析式上可能比后者更简单。例如Fisher信息矩阵可能是对角阵或者块对角阵，二观测信息矩阵可能不是。其次作为副产物，这两种算法都提供了极大似然估计的协方差矩阵。

广义线性模型（GLMs）中的最大似然估计（MLE）

首先，GLMs中的对数似然函数具有这样的形式 $l=∑i=1nlogf(yi;θi,ϕ)=∑i=1n(yiθi−b(θi))ai(ϕ)+∑i=1nc(yi,ϕ)l=\sum_{i=1}^{n}logf\left ( y_{i};\theta _{i} ,\phi \right )=\sum_{i=1}^{n}\frac{\left ( y_{i}\theta _{i}-b\left ( \theta _{i} \right ) \right )}{a_{i}\left ( \phi \right )}+\sum_{i=1}^{n}c\left ( y_{i} ,\phi \right )$ 其中， $β=(β1,β2,...,βp)T\beta =\left ( \beta _{1} ,\beta _{2},..., \beta _{p}\right )^{T}$ ， $βj\beta_{j}$ 的得分函数为 $Uj=∂l∂βj=∑i=1n(yi−b′(θi))ai(ϕ)∂θi∂βj=∑i=1n(yi−μ)ai(ϕ)∂θi∂βj（2）U_{j}=\frac{\partial l}{\partial \beta _{j}}=\sum_{i=1}^{n}\frac{\left ( y_{i}-b^{'}\left ( \theta _{i} \right ) \right )}{a_{i}\left ( \phi \right )}\frac{\partial \theta _{i}}{\partial \beta _{j}}=\sum_{i=1}^{n}\frac{\left ( y_{i}-\mu \right )}{a_{i}\left ( \phi \right )}\frac{\partial \theta _{i}}{\partial \beta _{j}}（2）$
其中， $μi=E(Yi)=b′(θi)，Var(Yi)=b′′(θi)a(ϕ)\mu _{i}=E\left ( Y_{i} \right )=b^{'}\left ( \theta _{i} \right )，Var\left ( Y_{i} \right )=b^{''}\left ( \theta _{i} \right )a\left ( \phi \right )$ ，我们使用链式法则进行差异化 $∂θi∂βj=∂θi∂μi∂μi∂βj\frac{\partial \theta _{i}}{\partial \beta _{j}}=\frac{\partial \theta _{i}}{\partial \mu _{i}}\frac{\partial \mu _{i}}{\partial \beta _{j}}$ 因为 $∂θi∂μi=1∂μi∂θi=1b′′(θi)=ai(ϕ)b′′(θi)ai(ϕ)=ai(ϕ)Var(Yi)\frac{\partial \theta _{i}}{\partial \mu _{i}}=\frac{1}{\frac{\partial \mu _{i}}{\partial \theta _{i}}}=\frac{1}{b^{''}\left ( \theta _{i} \right )}=\frac{a_{i}\left ( \phi \right )}{b^{''}\left ( \theta _{i} \right )a_{i}\left ( \phi \right )}=\frac{a_{i}\left ( \phi \right )}{Var\left ( Y_{i} \right )}$
并且 $∂μi∂βj=∂μi∂ηi∂ηi∂βj=∂μi∂ηixij\frac{\partial \mu _{i}}{\partial \beta _{j}}=\frac{\partial \mu _{i}}{\partial \eta _{i}}\frac{\partial \eta _{i}}{\partial \beta _{j}}=\frac{\partial \mu _{i}}{\partial \eta _{i}}x_{ij}$ 其中 $x_{ij}$ 是 $x_{i}的第j个分量$ ，我们知道 $∂θi∂βj=ai(ϕ)Var(Yi)∂μi∂ηjxij\frac{\partial \theta _{i}}{\partial \beta _{j}}=\frac{a_{i}\left ( \phi \right )}{Var\left ( Y_{i} \right )}\frac{\partial \mu _{i}}{\partial \eta _{j}}x_{ij}$ 因此（2）式就化为了 $Uj=∑i=1n[(yi−μi)Var(Yi)xij(∂μi∂ηi)]=∑i=1n(yi−μi)g′(μi)Vixij（3）U_{j}=\sum_{i=1}^{n}\left [ \frac{\left ( y_{i}-\mu _{i} \right )}{Var\left ( Y_{i} \right )}x_{ij} \left ( \frac{\partial \mu _{i}}{\partial \eta _{i}} \right )\right ]=\sum_{i=1}^{n}\frac{\left ( y_{i} -\mu _{i}\right )}{g^{'}\left ( \mu _{i} \right )V_{i}}x_{ij}（3）$ 其中 $Vi=Var(Yi)V_{i}=Var\left ( Y_{i} \right )$ ，并且 $∂μi∂ηi=1∂ηi∂μi=1g′(μi)\frac{\partial \mu _{i}}{\partial \eta _{i}}=\frac{1}{\frac{\partial \eta _{i}}{\partial \mu _{i}}}=\frac{1}{g^{'}\left ( \mu _{i}\right )}$ 由于 $ηi=g(μi)\eta _{i}=g\left ( \mu _{i} \right )$ ，因此 $β\beta$ 的得分向量是 $U≡U(β)=∑i=1n(yi−μi)g′(μi)Vixi（4）U\equiv U\left ( \beta \right )=\sum_{i=1}^{n}\frac{\left ( y_{i}-\mu _{i} \right )}{g^{'}\left ( \mu _{i} \right )V_{i}}x_{i}（4）$ 另一方面，（3）式对 $βj\beta_{j}$ 求偏导得 $∂2l∂βj∂βk=∂Uj∂βk=∑i=1n(−∂μi∂βk)1g′(μi)Vixij+∑i=1n(yi−μi)∂[1g′(μi)Vi]∂βkxij（5）\frac{\partial ^{2}l}{\partial \beta _{j}\partial \beta _{k}}=\frac{\partial U_{j}}{\partial \beta _{k}}=\sum_{i=1}^{n}\left ( -\frac{\partial \mu _{i}}{\partial \beta _{k}} \right )\frac{1}{g^{'}\left ( \mu _{i} \right )V_{i}}x_{ij}+\sum_{i=1}^{n}\left ( y_{i} -\mu _{i}\right )\frac{\partial \left [ \frac{1}{g^{'}\left ( \mu _{i} \right )V_{i}} \right ]}{\partial \beta _{k}}x_{ij}（5）$ 由于 $E(Yi−μi)=0E\left ( Y_{i} -\mu _{i}\right )=0$ ，所以（5）式的第二项在进行期望时就消失了。即Fisher信息阵的矩阵形式就变成了 $I(β)=E(∂2l∂β∂βT)=∑i=1n1g′(μi)2VixijxikI\left ( \beta \right )=E\left ( \frac{\partial ^{2}l}{\partial \beta \partial \beta ^{T}} \right )=\sum_{i=1}^{n}\frac{1}{g^{'} \left ( \mu _{i} \right )^{2}V_{i}}x_{ij}x_{ik}$ 因此当我们表示 $Wi=1g′(μi)2ViW_{i}=\frac{1}{g^{'}\left ( \mu _{i} \right )^{2}V_{i}}$ ，并且 $W=diag(W1,W2,...,Wn)=(W10⋯00W2⋯0⋮0⋱00⋯0Wn)W=diag\left ( W_{1},W_{2},...,W_{n} \right )=\begin{pmatrix} W_{1} &0 & \cdots & 0\\ 0& W_{2} & \cdots & 0\\ \vdots & 0 & \ddots & 0\\ 0& \cdots &0 & W_{n} \end{pmatrix}$ 则Fisher信息阵就可以表示为 $I(β)=XTWXI\left ( \beta \right )=X^{T}WX$ 令 $D=diag(g′(μ1),g′(μ2),...,g′(μn))D=diag\left ( g^{'}\left ( \mu _{1} \right ) ,g^{'}\left ( \mu _{2} \right ),...,g^{'}\left ( \mu _{n} \right )\right )$ ，这样（4）式就可以写成 $U=U(β)=XTWD(y−μ)U=U\left ( \beta \right )=X^{T}WD\left ( y-\mu \right )$

计算最大似然估计（MLE）参数 $β\beta$ 的算法

假设我们有一个估计 $β(m−1)\beta ^{\left ( m-1 \right )}$ ，基于这个估计我们计算 $μ(m−1)=μ(β(m−1))，W(m−1)=W(β(m−1))\mu ^{\left ( m-1 \right )}=\mu \left ( \beta ^{\left ( m-1 \right )} \right )，W^{\left ( m-1\right )}=W\left ( \beta ^{\left ( m-1 \right )} \right )$
并且有 $D(m−1)=D(β(m−1))D^{\left ( m-1 \right )}=D\left ( \beta ^{\left ( m-1 \right )} \right )$ 那么Fisher得分算法就会显示 $β\beta$ 的下一次迭代 $β(m)=β(m−1)+[I(β(m−1))]−1[U(β(m−1))]=β(m−1)+[XTW(M−1)X]−1[XTW(M−1)D(M−1)(y−μ(m−1))]\beta ^{\left ( m \right )}=\beta ^{\left ( m-1 \right )}+\left [ I\left ( \beta ^{\left ( m-1 \right )} \right ) \right ]^{-1}\left [ U\left ( \beta ^{\left ( m-1 \right )} \right ) \right ]=\beta ^{\left ( m-1 \right )}+\left [ X^{T} W^{\left ( M-1 \right )}X\right ]^{-1}\left [ X^{T} W^{\left ( M-1 \right )}D^{\left ( M-1 \right )}\left ( y-\mu ^{\left ( m-1 \right )} \right )\right ]$ 可以写成 $β(m)=[XTWX]−1XTW(m−1)[Xβ(m−1)+D(m−1)(y−μ(m−1))]\beta ^{\left ( m \right )}=\left [ X^{T}WX \right ]^{-1}X^{T}W^{\left ( m-1 \right )}\left [ X\beta ^{\left ( m-1 \right )}+D^{\left ( m-1 \right )} \left ( y-\mu ^{\left ( m-1 \right )} \right )\right ]$ 令 $Z(m−1)=Xβ(m−1)+D(m−1)(y−μ(m−1))Z^{\left ( m-1 \right )}=X\beta ^{\left ( m-1 \right )}+D^{\left ( m-1 \right )}\left ( y-\mu ^{\left ( m-1 \right )} \right )$ 然后它又可以被写成 $β(m)=(X(T)W(m−1)X)−1XTWm−1Z(m−1)（6）\beta ^{\left ( m \right )}=\left ( X^{\left ( T \right )}W^{\left ( m-1 \right )}X \right )^{-1}X^{T}W^{m-1}Z^{\left ( m-1 \right )}（6）$
注释（6）式意味着，给定参数 $β\beta$ 的解，我们需要计算“工作权重矩阵” $W$ 和“工作响应向量” $Z$ ，然后利用广义加权最小二乘法得到 $β\beta$ 的更新解。

广义线性模型实例解析

下表中的ARPI事物数据在协变量X的不同处观察到Y,并且数据是服从Poisson分布的。我们利用GLM来解决这个问题。

$Y_{i}$	2	3	6	7	8	9	10	12	15
$x_{i}$	-1	-1	0	0	0	0	1	1	1

数据即探索Y和X之间的关系。设 $Y_{i}$ 为变量 $y$ 的第 $i$ 个数，表示 $E(Yi)=μiE\left ( Y_{i} \right )=\mu _{i}$ 。我们通过建立关系 $g(μi)=xi′βg\left ( \mu _{i} \right )=x_{i}^{'}\beta$ 对于这个Poisson数据集，点则连接是对数连接函数。
$logμi=β0+β1xi=(1,xi)(β0β1)=xiTβlog\mu _{i}=\beta _{0}+\beta _{1}x_{i}=\left ( 1,x_{i} \right )\begin{pmatrix} \beta _{0}\\ \beta _{1}\end{pmatrix}=x_{i}^{T}\beta$ 接下来我们要来求 $W 和 Z$ 的表达式。
我们已知的条件有 $g′(μi)=1μig^{'}\left ( \mu _{i} \right )=\frac{1}{\mu _{i}}$ ，对于Poisson分布显然有 $Vi=E(Yi)=μiV_{i}=E\left ( Y_{i} \right )=\mu _{i}$ ，所以 $Wi=[(g′(μi)2)Vi]−1=exp{xiTβ}W_{i}=\left [ \left ( g^{'}\left ( \mu _{i} \right ) ^{2}\right )V_{i} \right ]^{-1}=exp\left \{ x_{i}^{T} \beta \right \}$ 并且 $Zi=xiTβ+g′(μi)(yi−μi)=xiTβ+(yi−μi)μiZ_{i}=x_{i}^{T}\beta +g^{'}\left ( \mu _{i} \right )\left ( y_{i} -\mu _{i}\right )=x_{i}^{T}\beta +\frac{\left ( y_{i}-\mu _{i} \right )}{\mu _{i}}$
我们选择 $β的初始值β0=2，β1=1\beta的初始值\beta_{0}=2，\beta_{1}=1$ 。结合Fisher迭代算法，代入数据。这个过程一直持续到收敛。结果如下表

m	0	1	2	3	4
$β0m\beta _{0}^{m}$	2	1.9150	1.8902	1.8892	1.8892
$β1m\beta _{1}^{m}$	1	0.7235	1.8902	1.8892	1.8892

因此 $β的MLE是β0=1.8892，β1=0.6697\beta的MLE是\beta _{0}=1.8892，\beta _{1}=0.6697$

R语言代码

y <- c(2,3,6,7,8,9,10,12,15); 
x <- c(-1,-1,0,0,0,0,1,1,1)
X <- cbind(rep(1,9),x); beta_0 <- c(2,1)
for (i in 1:100){
beta <- beta_0
eta <- X %*% beta
mu <- exp(eta)
W <- diag(as.vector(mu))
Z <- X %*% beta + ((y-mu)*mu^(-1))
XWX <- t(X) %*% W %*% X
XWZ <- t(X) %*% W %*% Z
Cov <- solve(XWX)
beta_0 <- Cov %*% XWZ}
testdata<-data.frame(y,x)
summary(glm(y~x,family=poisson,data=testdata))

2048 AI社区

有“AI”的1024 = 2048，欢迎大家加入2048 AI社区

更多推荐

SSH 命令从新手入门到安全运维

2048 AI社区

CAD-Assistant 闭环逻辑详解：无训练也能精准建模的核心密码

CAD-Assistant 的闭环逻辑，本质是用 “大模型的推理能力 + 工具的精确性 + 环境的反馈机制”，重构了 CAD 设计的自动化流程。其核心价值在于：零训练成本：无需标注 CAD 数据集，仅通过工具调用和反馈就能适配多样场景；高精度落地：闭环迭代确保结果与需求一致，解决了 VLLM “光说不练” 的痛点；低使用门槛：设计师无需编写代码，仅通过草图、文本就能生成专业 CAD 模型。