泊松回归代码python 泊松回归适用条件

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

泊松回归

定义

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

• e的概念
  $\lim_{n\to\infty}(1-\frac{1}{n})^n=e \\ \lim_{n\to\infty}(1-\frac{\lambda}{n})^n=e^\lambda \\$
• 二项式分布
  $P(X=k)=\binom{n}{k}p^k(1-p)^{n-k}$
  如果我们令$p=\lambda/n$,当$n\to\infty$时$P$的极限是:
  $\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\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{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}$
  再做指数变换:
  $e^{log(Y;\lambda)}=f(Y;\lambda)=\exp\{klog(\lambda)-\lambda-log(k!) \}$
  其中$log(Y)$称为链接函数,由此得到了泊松回归模型的形式.那么为什么说泊松回归模型是一种广义线性模型(GLM)呢？首先考虑线性回归模型:
  $y=X\beta+e \\ \eta=X\beta\\ \mu=E(y)$
  在广义线性模型中,我们不再要求$y$服从$N(\mu,\sigma^2)$而是服从于指数分布族,例如Normal;Poisson;Gamma;Binomial.现证明泊松回归属于指数分布族,首先考虑指数分布族的定义：
  $f(y;\theta,\phi)=\exp\left\{ \frac{(y\theta-b(\theta))}{a(\phi)} +c(y,\phi)\right\}$
  现在我们回顾一下泊松回归模型的形式:
  $f(Y;\lambda)=\exp\{klog(\lambda)-\lambda-log(k!) \}$
  记$\theta=log(\lambda),\lambda=\exp\{\theta\}$,则有：
  $\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}$
  所以泊松回归模型也属于指数分布族

参数估计

• 似然函数
  再给出了模型的形式以后,进一步需要对参数进行估计,采用极大似然估计法,从总体$(Y,X)$中抽取容量为$n$的随机样本,此时有似然函数:
  $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(\beta)=\sum_{i=1}^n[y_ix_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_ix_i-x_i\exp(x_i\beta_i)]=0$
  上式方程组一般情况下并没有解析解,但我们可以通过牛顿拉夫逊迭代法求解:
  $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}$
  其中$\beta_0=1,x_0=1$.则$F(\beta_i)$关于$\beta$的雅克比矩阵:
  $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.$
  此时有Newton-Raphson算法：
  $\beta^{(m+1)}=\beta^{(m)}-[J(\beta^{(m)})]^{-1}F(\beta^{(m)}),m=0,1,\cdots$
  在这个迭代过程中,需要给定$\beta$的初值和精度$\varepsilon$,不断计算上述过程直至$|\beta^{(m+1)}-\beta^{(m)}|<\varepsilon$收敛后结束.
  同时附上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);