最近在研究GWPR,参考了很多广义线性模型,特别是泊松回归的相关内容,知识琐碎且繁杂,做个笔记。

泊松回归

定义

泊松回归(Poisson regression)是用来为计数资料和列联表建模的一种回归分析.泊松回归假设反应变量Y是泊松分布,并假设它期望值的对数可被未知参数的线性组合建模.泊松回归模型有时(特别是当用作列联表模型时)又被称作对数-线性模型.需要注意的是,对数线性模型和泊松回归模型并不完全相同,通常对数线性回归的响应变量是连续的,而泊松回归则是离散的.再给出泊松回归模型的形式之前,我们先考虑几个概念:

e的概念
lim⁡n→∞(1−1n)n=elim⁡n→∞(1−λn)n=eλ\lim_{n\to\infty}(1-\frac{1}{n})^n=e \\ \lim_{n\to\infty}(1-\frac{\lambda}{n})^n=e^\lambda \\ n→∞lim(1−n1)n=en→∞lim(1−nλ)n=eλ
二项式分布
P(X=k)=(nk)pk(1−p)n−kP(X=k)=\binom{n}{k}p^k(1-p)^{n-k} P(X=k)=(kn)pk(1−p)n−k
如果我们令p=λ/np=\lambda/np=λ/n,当n→∞n\to\inftyn→∞时PPP的极限是:
lim⁡n→∞P(X=k)=lim⁡n→∞(nk)pk(1−p)n−k=lim⁡n→∞n!(n−k)!k!(λn)k(1−λn)n−k=lim⁡n→∞[n!nk(n−k)!](λkk!)(1−λn)n⏟→exp⁡(−λ)(1−λn)−k⏟→1=lim⁡n→∞[(1−1n)(1−2n)⋯(1−k−1n)](λkk!)(1−1n)n(1−λn)−k=(λkk!)exp⁡{−λ}\begin{aligned} \lim_{n\to\infty}P(X=k)&=\lim_{n\to\infty}\binom{n}{k}p^k(1-p)^{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}\left[\frac{n!}{n^k(n-k)!}\right] \left( \frac{\lambda^k}{k!}\right) \underbrace{ \left( 1-\frac{\lambda}{n} \right) ^n }_{\to\exp(-\lambda)}\underbrace{\left( 1-\frac{\lambda}{n} \right) ^{-k} }_{\to1}\\ &=\lim_{n\to\infty} \left[ \left( 1-\frac{1}{n} \right) \left( 1-\frac{2}{n} \right)\cdots \left( 1-\frac{k-1}{n} \right)\right]\left( \frac{\lambda^k}{k!} \right) \left( 1-\frac{1}{n} \right)^n \left( 1-\frac{\lambda}{n} \right)^{-k} \\ &=\left( \frac{\lambda^k}{k!} \right)\exp\{-\lambda\} \end{aligned} n→∞limP(X=k)=n→∞lim(kn)pk(1−p)n−k=n→∞lim(n−k)!k!n!(nλ)k(1−nλ)n−k=n→∞lim[nk(n−k)!n!](k!λk)→exp(−λ)(1−nλ)n→1(1−nλ)−k=n→∞lim[(1−n1)(1−n2)⋯(1−nk−1)](k!λk)(1−n1)n(1−nλ)−k=(k!λk)exp{−λ}
这也说明当n→∞n\to\inftyn→∞时,二项式分布可以近似至泊松分布
泊松分布与泊松回归
P(X=k)=λkk!e−λ,k=0,1,⋯P(X=k)=\frac{\lambda^k}{k!}e^{-\lambda},k=0,1,\cdots P(X=k)=k!λke−λ,k=0,1,⋯
其中有E(X)=λ,Var(X)=λE(X)=\lambda,Var(X)=\lambdaE(X)=λ,Var(X)=λ,为导出泊松回归的形式,我们记上式为YYY,则有:
Y=λkk!e−λ,k=0,1,⋯Y=\frac{\lambda^k}{k!}e^{-\lambda},k=0,1,\cdots Y=k!λke−λ,k=0,1,⋯
做对数变换:
log(Y;λ)=log(e−λλk)−log(k!)=log(e−λ)+log(λk)−log(k!)=−λ+klog(λ)−log(k!)\begin{aligned} log(Y;\lambda)&=log(e^{-\lambda}\lambda^k)-log(k!)\\ &=log(e^{-\lambda})+log(\lambda^k)-log(k!)\\ &=-\lambda+klog(\lambda)-log(k!) \\ \end{aligned} log(Y;λ)=log(e−λλk)−log(k!)=log(e−λ)+log(λk)−log(k!)=−λ+klog(λ)−log(k!)
再做指数变换:
elog(Y;λ)=f(Y;λ)=exp⁡{klog(λ)−λ−log(k!)}e^{log(Y;\lambda)}=f(Y;\lambda)=\exp\{klog(\lambda)-\lambda-log(k!) \} elog(Y;λ)=f(Y;λ)=exp{klog(λ)−λ−log(k!)}
其中log(Y)log(Y)log(Y)称为链接函数,由此得到了泊松回归模型的形式.那么为什么说泊松回归模型是一种广义线性模型(GLM)呢？首先考虑线性回归模型:
y=Xβ+eη=Xβμ=E(y)y=X\beta+e \\ \eta=X\beta\\ \mu=E(y) y=Xβ+eη=Xβμ=E(y)
在广义线性模型中,我们不再要求yyy服从N(μ,σ2)N(\mu,\sigma^2)N(μ,σ2)而是服从于指数分布族,例如Normal;Poisson;Gamma;Binomial.现证明泊松回归属于指数分布族,首先考虑指数分布族的定义：
f(y;θ,ϕ)=exp⁡{(yθ−b(θ))a(ϕ)+c(y,ϕ)}f(y;\theta,\phi)=\exp\left\{ \frac{(y\theta-b(\theta))}{a(\phi)} +c(y,\phi)\right\} f(y;θ,ϕ)=exp{a(ϕ)(yθ−b(θ))+c(y,ϕ)}
现在我们回顾一下泊松回归模型的形式:
f(Y;λ)=exp⁡{klog(λ)−λ−log(k!)}f(Y;\lambda)=\exp\{klog(\lambda)-\lambda-log(k!) \} f(Y;λ)=exp{klog(λ)−λ−log(k!)}
记θ=log(λ),λ=exp⁡{θ}\theta=log(\lambda),\lambda=\exp\{\theta\}θ=log(λ),λ=exp{θ},则有：
f(Y=k;θ,ϕ)=exp⁡{yθ−exp⁡{θ}−log(k!)}=exp⁡{yθ−exp⁡{θ}1+(−)log(k!)}\begin{aligned} f(Y=k;\theta,\phi)&=\exp\{y\theta-\exp\{\theta\}-log(k!) \} \\ &=\exp\left\{ \frac{y\theta-\exp\{\theta\}}{1}+(-)log(k!)\right\} \end{aligned} f(Y=k;θ,ϕ)=exp{yθ−exp{θ}−log(k!)}=exp{1yθ−exp{θ}+(−)log(k!)}
所以泊松回归模型也属于指数分布族

参数估计

似然函数
再给出了模型的形式以后,进一步需要对参数进行估计,采用极大似然估计法,从总体(Y,X)(Y,X)(Y,X)中抽取容量为nnn的随机样本,此时有似然函数:
L(β)=∏i=1nf(β;yi,xi)=∏i=1nexp⁡{yiXiTβi−exp⁡(XiTβi)−log⁡(yi!)}L(\beta)=\prod_{i=1}^nf(\beta;y_i,x_i)=\prod_{i=1}^n\exp\{ y_iX_i^{\rm{T}}\beta_i-\exp(X_i^{\rm{T}}\beta_i)-\log(y_i!)\} L(β)=i=1∏nf(β;yi,xi)=i=1∏nexp{yiXiTβi−exp(XiTβi)−log(yi!)}
对数似然:
l(β)=∑i=1n[yixiβi−exp⁡(xiβi)−log⁡(yi!)]l(\beta)=\sum_{i=1}^n[y_ix_i\beta_i-\exp(x_i\beta_i)-\log(y_i!)] l(β)=i=1∑n[yixiβi−exp(xiβi)−log(yi!)]
对βi\beta_iβi求偏导:
∂l(β)∂βi=∑i=1n[yixi−xiexp⁡(xiβi)]=0\frac{\partial l(\beta)}{\partial\beta_i}=\sum_{i=1}^n[y_ix_i-x_i\exp(x_i\beta_i)]=0 ∂βi∂l(β)=i=1∑n[yixi−xiexp(xiβi)]=0
上式方程组一般情况下并没有解析解,但我们可以通过牛顿拉夫逊迭代法求解:
F(β)=[∑i=1n[yi−exp⁡(xiβi)]∑i=1n[yixi1−xi1exp⁡(xiβi)]⋮∑i=1n[yixiq−xiqexp⁡(xiβi)]]F(\beta)= \begin{bmatrix} \sum_{i=1}^n [y_i-\exp(x_i\beta_i)] \\ \sum_{i=1}^n [y_ix_{i1}-x_{i1}\exp(x_i\beta_i)] \\ \vdots\\ \sum_{i=1}^n [y_ix_{iq}-x_{iq}\exp(x_i\beta_i)] \end{bmatrix} F(β)=⎣⎢⎢⎢⎡∑i=1n[yi−exp(xiβi)]∑i=1n[yixi1−xi1exp(xiβi)]⋮∑i=1n[yixiq−xiqexp(xiβi)]⎦⎥⎥⎥⎤
其中β0=1,x0=1\beta_0=1,x_0=1β0=1,x0=1.则F(βi)F(\beta_i)F(βi)关于β\betaβ的雅克比矩阵:
J(β)=∂2l(β)∂βk∂βj=−∑i=1nxixikexp⁡(XTβ),k=0,1,⋯,q;j=0,1,⋯,q.J(\beta)=\frac{\partial^2l(\beta)}{\partial\beta_k\partial\beta_j}=-\sum_{i=1}^nx_ix_{ik}\exp(X^T\beta),k=0,1,\cdots,q;j=0,1,\cdots,q. J(β)=∂βk∂βj∂2l(β)=−i=1∑nxixikexp(XTβ),k=0,1,⋯,q;j=0,1,⋯,q.
此时有Newton-Raphson算法：
β(m+1)=β(m)−[J(β(m))]−1F(β(m)),m=0,1,⋯\beta^{(m+1)}=\beta^{(m)}-[J(\beta^{(m)})]^{-1}F(\beta^{(m)}),m=0,1,\cdots β(m+1)=β(m)−[J(β(m))]−1F(β(m)),m=0,1,⋯
在这个迭代过程中,需要给定β\betaβ的初值和精度ε\varepsilonε,不断计算上述过程直至∣β(m+1)−β(m)∣<ε|\beta^{(m+1)}-\beta^{(m)}|<\varepsilon∣β(m+1)−β(m)∣<ε收敛后结束.
同时附上MATLAB代码

function F = PoissionRegressopt(b,Y,X)
n = length(Y);
F = 0;
for k = 1:nF = F + Y(k)*X(k,:)*b - exp(X(k,:)*b);% - factorial(Y(k));
end
F = - F;function F = PoissionF(b,Y,X)
n = length(Y);
F = zeros(size(b));
for k = 1:nF = F + Y(k)*X(k,:)'- exp(X(k,:)*b)*X(k,:)';
endfunction JM  = PoissionJM(b,Y,X)
n = length(Y);
JM  = zeros(size(b,1));
for k = 1:nJM = JM +  exp(X(k,:)*b)*X(k,:)'*X(k,:);
endfunction [ bm fv1,fv2] = PoissionNR(bm0,Y,X)
itermax = 30;
errstol = 1e-4;
iters = 0;
deltabm = ones(size(bm0));
bm1 = bm0 + deltabm;
while (iters<itermax)||(max(abs(deltabm))>errstol)deltabm  = pinv(PoissionJM(bm0,Y,X))*PoissionF(bm0,Y,X);bm1 = bm0 + deltabm;bm0  = bm1;   iters = iters +1;
end
bm = bm0;
fv1 = PoissionF(bm,Y,X);
fv2 = PoissionRegressopt(bm,Y,X);

广义线性模型之泊松回归相关推荐

Python使用sklearn构建广义线性模型：泊松回归（Poisson regression）实战
Python使用sklearn构建广义线性模型:泊松回归(Poisson regression)实战目录 Python使用sklearn构建广义线性模型:泊松回归(Poisson regressio ...
Python使用sklearn构建广义线性模型：gamma回归（Gamma regression）实战
Python使用sklearn构建广义线性模型:gamma回归(Gamma regression)实战目录 Python使用sklearn构建广义线性模型:gamma回归(Gamma regress ...
Python使用sklearn构建广义线性模型：Tweedie回归（Tweedie regression）实战
Python使用sklearn构建广义线性模型:Tweedie回归(Tweedie regression)实战目录 Python使用sklearn构建广义线性模型:Tweedie回归(Tweedie ...
广义线性模型（逻辑回归、泊松回归）
线性回归模型也并不适用于所有情况,有些结果可能包含而元数据(比如正面与反面)或者计数数据,广义线性模型可用于解释这类数据,使用的仍然是自变量的线性组合. 目录逻辑回归使用statsmodels 使 ...
stats | 广义线性模型（二）——泊松回归
本篇介绍一种常见的广义线性模型:泊松回归.泊松分布是离散型分布,它的概率分布函数如下: 写成指数族分布的形式如下: 对照指数族分布的通式: 可得, 广义线性模型假设与解释变量存在线性关系,即又因为泊 ...
对数线性模型之一(逻辑回归), 广义线性模型学习总结
经典线性模型自变量的线性预测就是因变量的估计值. 广义线性模型:自变量的线性预测的函数是因变量的估计值.常见的广义线性模型有:probit模型.poisson模型.对数线性模型等等.对数线性模型里有: ...
大数据分析R中泊松回归模型实例
如果您知道如何以及何时使用泊松回归,它可能是一个非常有用的工具.在大数据分析R中泊松回归模型实例中,我们将深入研究泊松回归,它是什么以及R程序员如何在现实世界中使用它. 具体来说,我们将介绍: 1)泊 ...
调整泊松回归中的过度分散
广义线性模型(Generalized Linear Models) The Poisson regression model naturally arises when we want to mode ...
泊松分布与泊松回归模型
泊松分布 Poisson分布(法语:loi de Poisson,英语:Poisson distribution,译名有泊松分布.普阿松分布.卜瓦松分布.布瓦松分布.布阿松分布.波以松分布.卜氏分配等 ...
R语言广义线性模型函数GLM、广义线性模型（Generalized linear models）、GLM函数的语法形式、glm模型常用函数、常用连接函数、逻辑回归、泊松回归、系数解读、过散度分析
R语言广义线性模型函数GLM.广义线性模型(Generalized linear models).GLM函数的语法形式.glm模型常用函数.常用连接函数.逻辑回归.泊松回归.系数解读.过散度分析目录

广义线性模型之泊松回归

泊松回归

定义

参数估计

广义线性模型之泊松回归相关推荐

最新文章

热门文章