大家好，欢迎来到IT知识分享网。

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

泊松回归

定义

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

• e的概念
  $$\begin{gathered} \underset{n\to \infty }{\lim }(1-\frac{1}{n}{)}^n=e \\ \underset{n\to \infty }{\lim }(1-\frac{\lambda }{n}{)}^n=e^{\lambda } \end{gathered}$$

• 二项式分布
  $$P(X=k)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^k(1-p{)}^{n-k}$$
  如果我们令$p=\lambda /n$,当$n\to \infty$时$P$的极限是:
  lim ⁡ n → ∞ P ( X = k ) = lim ⁡ n → ∞ ( n k ) p k ( 1 − p ) n − k = lim ⁡ n → ∞ n ! ( n − k ) ! k ! ( λ n ) k ( 1 − λ n ) n − k = lim ⁡ n → ∞ [ n ! n k ( n − k ) ! ] ( λ k k ! ) ( 1 − λ n ) n ⏟ → exp ⁡ ( − λ ) ( 1 − λ n ) − k ⏟ → 1 = lim ⁡ n → ∞ [ ( 1 − 1 n ) ( 1 − 2 n ) ⋯ ( 1 − k − 1 n ) ] ( λ k k ! ) ( 1 − 1 n ) n ( 1 − λ n ) − k = ( λ k k ! ) 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→∞lim​P(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\to \infty$时,二项式分布可以近似至泊松分布

• 泊松分布与泊松回归
  $$P(X=k)=\frac{{\lambda }^k}{k!}{e}^{-\lambda },k=0,1,\cdots$$
  其中有$E(X)=\lambda ,Var(X)=\lambda$,为导出泊松回归的形式,我们记上式为$Y$,则有:
  $$Y=\frac{{\lambda }^k}{k!}{e}^{-\lambda },k=0,1,\cdots$$
  做对数变换:
  $$\begin{array}{rl}log(Y;\lambda )& =log(e^{-\lambda }{\lambda }^k)-log(k!)\\ & =log(e^{-\lambda })+log({\lambda }^k)-log(k!)\\ & =-\lambda +klog(\lambda )-log(k!)\end{array}$$
  再做指数变换:
  e l o g ( Y ; λ ) = f ( Y ; λ ) = exp ⁡ { k l o g ( λ ) − λ − l o g ( k ! ) } e^{log(Y;\lambda)}=f(Y;\lambda)=\exp\{klog(\lambda)-\lambda-log(k!) \} elog(Y;λ)=f(Y;λ)=exp{
  klog(λ)−λ−log(k!)}
  其中$log(Y)$称为链接函数,由此得到了泊松回归模型的形式.那么为什么说泊松回归模型是一种广义线性模型(GLM)呢？首先考虑线性回归模型:
  $$\begin{gathered} y=X\beta +e \\ \eta =X\beta \\ \mu =E(y) \end{gathered}$$
  在广义线性模型中,我们不再要求$y$服从$N(\mu ,{\sigma }^2)$而是服从于指数分布族,例如Normal;Poisson;Gamma;Binomial.现证明泊松回归属于指数分布族,首先考虑指数分布族的定义：
  $$f(y;\theta ,\varphi )=\exp \left\{\frac{(y\theta -b(\theta ))}{a(\varphi )}+c(y,\varphi )\right\}$$
  现在我们回顾一下泊松回归模型的形式:
  f ( Y ; λ ) = exp ⁡ { k l o g ( λ ) − λ − l o g ( k ! ) } f(Y;\lambda)=\exp\{klog(\lambda)-\lambda-log(k!) \} f(Y;λ)=exp{
  klog(λ)−λ−log(k!)}
  记 θ = l o g ( λ ) , λ = exp ⁡ { θ } \theta=log(\lambda),\lambda=\exp\{\theta\} θ=log(λ),λ=exp{
  θ},则有：
  f ( Y = k ; θ , ϕ ) = exp ⁡ { y θ − exp ⁡ { θ } − l o g ( k ! ) } = exp ⁡ { y θ − exp ⁡ { θ } 1 + ( − ) l o g ( 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)$中抽取容量为$n$的随机样本,此时有似然函数:
  L ( β ) = ∏ i = 1 n f ( β ; y i , x i ) = ∏ i = 1 n exp ⁡ { y i X i T β i − exp ⁡ ( X i T β i ) − log ⁡ ( y i ! ) } 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∏n​f(β;yi​,xi​)=i=1∏n​exp{
  yi​XiT​βi​−exp(XiT​βi​)−log(yi​!)}
  对数似然:
  $$l(\beta )=\sum _{i=1}^n[y_i{x}_i{\beta }_i-\exp (x_i{\beta }_i)-\log (y_i!)]$$
  对${\beta }_i$求偏导:
  $$\frac{\partial l(\beta )}{\partial {\beta }_i}=\sum _{i=1}^n[y_i{x}_i-x_i\exp (x_i{\beta }_i)]=0$$
  上式方程组一般情况下并没有解析解,但我们可以通过牛顿拉夫逊迭代法求解:
  $$F(\beta )=\left[\begin{array}{c}{\sum }_{i=1}^n[y_i-\exp (x_i{\beta }_i)]\\ {\sum }_{i=1}^n[y_i{x}_{i1}-x_{i1}\exp (x_i{\beta }_i)]\\ ⋮\\ {\sum }_{i=1}^n[y_i{x}_{iq}-x_{iq}\exp (x_i{\beta }_i)]\end{array}\right]$$
  其中${\beta }_0=1,x_0=1$.则$F({\beta }_i)$关于$\beta$的雅克比矩阵:
  $$J(\beta )=\frac{{\partial }^{2}l(\beta )}{\partial {\beta }_k\partial {\beta }_j}=-\sum _{i=1}^n{x}_i{x}_{ik}\exp (X^T\beta ),k=0,1,\cdots ,q;j=0,1,\cdots ,q.$$
  此时有Newton-Raphson算法：
  $${\beta }^{(m+1)}={\beta }^{(m)}-[J({\beta }^{(m)}){]}^{-1}F({\beta }^{(m)}),m=0,1,\cdots$$
  在这个迭代过程中,需要给定$\beta$的初值和精度$\epsilon$,不断计算上述过程直至$|{\beta }^{(m+1)}-{\beta }^{(m)}|<\epsilon$收敛后结束.
  同时附上MATLAB代码

function F = PoissionRegressopt(b,Y,X)
n = length(Y);
F = 0;
for k = 1:n
    F = 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:n
    F = F + Y(k)*X(k,:)'- exp(X(k,:)*b)*X(k,:)';
end


function JM  = PoissionJM(b,Y,X)
n = length(Y);
JM  = zeros(size(b,1));
for k = 1:n
    JM = JM +  exp(X(k,:)*b)*X(k,:)'*X(k,:);
end

function [ 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);