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

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

广义线性模型（GLM）定义

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

连接函数介绍

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

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

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

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

假设我们对未知参数θ\thetaθ的对数似然函数，比如说l(θ)l\left ( \theta \right )l(θ)我们想找出θ\thetaθ的最大似然估计（MLE） θ^\hat{\theta }θ^，即θ^≡argmaxθ⊂Ω{l(θ)}\hat{\theta }\equiv arg \underset{\theta \subset \Omega }{max}\left \{ l\left ( \theta \right ) \right \}θ^≡argθ⊂Ωmax{l(θ)}
这是估计方程的解∂l(θ)∂θ=0\frac{\partial l\left ( \theta \right )}{\partial \theta }=0∂θ∂l(θ)=0
1、在这种情况下，例如，对于正态分布参数θ\thetaθ的最大似然估计 θ^\hat{\theta }θ^可以有一个明确的数学表达式（μ=1n∑i=1nlnxi）（\mu =\frac{1}{n}\sum_{i=1}^{n}lnx_{i}）（μ=n1∑i=1nlnxi）
2、一般来说，最大似然估计θ\thetaθ没有没有明确的数学解。相反，需要某些数值优化方法。
3、统计学中最常用的两种优化方法是Newton-Raphson算法和Fisher得分算法，他们都涉及计算l(θ)l\left ( \theta \right )l(θ)对θ\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=θm−1+[−l′′(θ(m−1))]−1[l′(θ(m−1))]（1）
其中m=1,2,...m=1,2,...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))=∂θ∂l(θ)∣θ=θ(m−1)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 )}}l′′(θ(m−1))=∂θ∂θT∂2l(θ)∣θ=θ(m−1)是p×1p\times 1p×1和p×pp\times pp×p的向量和矩阵（p是θ的维数）（p是\theta的维数）（p是θ的维数）
注1 l′(θ)被称为θ的得分函数。−l′′(θ)被称为θ的观测信息矩阵l^{'}\left ( \theta \right )被称为\theta的得分函数。-l^{''}\left ( \theta \right )被称为\theta的观测信息矩阵l′(θ)被称为θ的得分函数。−l′′(θ)被称为θ的观测信息矩阵
注2 算法（1）需要初始值，例如θ0\theta ^{0}θ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−1)∥∥∥∥θ(m)−θ(m−1)∥∥≤10−5则迭代停止。θ(m)\theta ^{\left ( m \right )}θ(m)可以认为是最大似然估计 θ^\hat{\theta }θ^。

Fisher得分算法

Fisher得分算法与Newton-Raphson算法相同，只是（1）式中的观测矩阵−l′′(θ)-l^{''}\left ( \theta \right )−l′′(θ)被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 )dYI(θ)=E[−l′′(θ)]=−∫l′′(θ∣Y)fY(Y∣θ)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 )l=i=1∑nlogf(yi;θi,ϕ)=i=1∑nai(ϕ)(yiθi−b(θi))+i=1∑nc(yi,ϕ)其中，β=(β1,β2,...,βp)T\beta =\left ( \beta _{1} ,\beta _{2},..., \beta _{p}\right )^{T}β=(β1,β2,...,βp)T，βj\beta_{j}β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）Uj=∂βj∂l=i=1∑nai(ϕ)(yi−b′(θi))∂βj∂θi=i=1∑nai(ϕ)(yi−μ)∂βj∂θi（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=E(Yi)=b′(θi)，Var(Yi)=b′′(θi)a(ϕ)，我们使用链式法则进行差异化∂θ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}}∂βj∂θi=∂μi∂θi∂βj∂μi因为∂θ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∂θi=∂θi∂μi1=b′′(θi)1=b′′(θi)ai(ϕ)ai(ϕ)=Var(Yi)ai(ϕ)
并且∂μ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}∂βj∂μi=∂ηi∂μi∂βj∂ηi=∂ηi∂μixij其中xijx_{ij}xij是xi的第j个分量x_{i}的第j个分量xi的第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}∂βj∂θi=Var(Yi)ai(ϕ)∂ηj∂μixij因此（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）Uj=i=1∑n[Var(Yi)(yi−μi)xij(∂ηi∂μi)]=i=1∑ng′(μi)Vi(yi−μi)xij（3）其中Vi=Var(Yi)V_{i}=Var\left ( Y_{i} \right )Vi=Var(Yi)，并且∂μ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∂μi=∂μi∂ηi1=g′(μi)1由于ηi=g(μi)\eta _{i}=g\left ( \mu _{i} \right )ηi=g(μi)，因此β\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）U≡U(β)=i=1∑ng′(μi)Vi(yi−μi)xi（4）另一方面，（3）式对βj\beta_{j}β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）∂βj∂βk∂2l=∂βk∂Uj=i=1∑n(−∂βk∂μi)g′(μi)Vi1xij+i=1∑n(yi−μi)∂βk∂[g′(μi)Vi1]xij（5）由于E(Yi−μi)=0E\left ( Y_{i} -\mu _{i}\right )=0E(Yi−μi)=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}I(β)=E(∂β∂βT∂2l)=i=1∑ng′(μi)2Vi1xijxik因此当我们表示Wi=1g′(μi)2ViW_{i}=\frac{1}{g^{'}\left ( \mu _{i} \right )^{2}V_{i}}Wi=g′(μi)2Vi1，并且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}W=diag(W1,W2,...,Wn)=⎝⎜⎜⎜⎛W10⋮00W20⋯⋯⋯⋱0000Wn⎠⎟⎟⎟⎞则Fisher信息阵就可以表示为I(β)=XTWXI\left ( \beta \right )=X^{T}WXI(β)=XTWX令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 )D=diag(g′(μ1),g′(μ2),...,g′(μn))，这样（4）式就可以写成U=U(β)=XTWD(y−μ)U=U\left ( \beta \right )=X^{T}WD\left ( y-\mu \right )U=U(β)=XTWD(y−μ)

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

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

广义线性模型实例解析

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

YiY_{i}Yi	2	3	6	7	8	9	10	12	15
xix_{i}xi	-1	-1	0	0	0	0	1	1	1

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

m	0	1	2	3	4
β0m\beta _{0}^{m}β0m	2	1.9150	1.8902	1.8892	1.8892
β1m\beta _{1}^{m}β1m	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β的MLE是β0=1.8892，β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))