Poission 回归参数最大似然估计的计算

Poisson 回归参数最大似然估计的计算
1 Possion 回归模型的定义
假设因变量Y 是一个服从Poission 分布的随机变量，12,,,q x x x 是影响Y 的k 个因素，12[1,,,,]T q X x x x =是协变量向量，01[,,,]T q βββ=β是回归参数向量，则Y 关于x 的k 元Poission 回归模型定义为
()
{}exp(()),0,1,2,.!k X P Y k X X k k λλ==-= (1)
其中()exp()0.T X X λ=>β
2 参数估计
我们用最大似然估计方法去求模型的参数。

假设从总体(,)Y X 中抽取一个容量为n 的随机样本1122(,),(,),
,(,)n n y X y X y X ，其中12[1,,,,],1,2,
,T k k k kq X x x x k n ==，则有似然函数为 11
exp ()(){}exp(exp())!k y T n n T k k k k k k k X L P Y y X X y =====-∏∏βββ (2) 两边取对数，整理可得
1ln ()exp()ln(!)n
T T k k k k k L y X X y =⎡⎤=--⎣⎦∑βββ (3)
为研究方便，以下不妨记01k x =。

为求式(3)的最大值点，即最大似然估计，可求对数似然函数ln ()L β关于β的似然方程组为
1
l n ()[e x p ()]n T k k i k i k k i L y x x X β=∂=-∂∑ββ，0,1,,.i q = (4)
具体形式为 1011111ln ()[exp()]0ln ()[exp()]0ln ()[exp()]0n T k k k n T k k k k k n T k kq kq k k k L y X L y x x X L y x x X βββ===∂⎧=-=⎪∂⎪⎪∂=-=⎪
∂⎨⎪⎪⎪∂=-=⎪∂⎩
∑∑∑ββββββ (5) 式(5) 为非线性方程组，一般情况下没有解析解，可以用Newton-Raphson 迭代方法求其数值解，令
11111[exp()][exp()]()[exp()]n T k k k n T k k k k k n T k kq kq k k y X y x x X F y x x X ===⎡⎤-⎢⎥⎢⎥⎢⎥-⎢⎥=⎢⎥⎢⎥⎢⎥
⎢⎥-⎢⎥⎣⎦
∑∑∑ββββ (6) 则()F β关于β的Jacobian 矩阵为
21ln ()()exp(),0,1,,,0,1,,.n T ki kj k i j L J x x X i q j q βββ=∂==-==∂∂∑ββ (7)
具体形式为
1111211111121111exp()exp()exp()exp()exp()exp()()exp()exp()exp()n n n T T T k k k kq k k k k n n n T T T k k k k k kq k k k k n n n T T T kq k kq k k kq k k k k X x X x X x X x X x x X J x X x x X x X =========⎡⎤---⎢⎥⎢⎥⎢⎥---⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥---⎢⎣⎦∑∑∑∑∑∑∑∑∑ββββββββββ⎥ (7)
对应的向量形式为
1()exp()n
T T k k k
k J X X X ==-∑ββ (7’) 根据Newton-Raphson 方法的原理，可得参数β迭代公式为
1(1)()()()()(),0,1,2,.m m m m J F m -+⎡⎤=-=⎣⎦ββββ (8)
算法如下：
Step 1: 给定参数β的初值参数(0)β和误差容许精度ε，令0m =；
Step 2：计算1(1)()()()()(),0,1,2,.m m m m J F m -+⎡⎤=-=⎣⎦ββββ；
Step 3: 若()()m F ε<β，即满足容许的精度，则结束，否则更新参数()(1)m m +=ββ，
1m m =+，转至Step2.
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);
附录1:
>> b =glmfit(X0,Y,'poisson', 'log')
b =
1.5043
0.4518
0.3578
0.2388
可以看到，结果一致。

比文献【1】中的结果要好一点
参考文献
【1】茆诗松主编. 统计手册[M]. 北京: 科学出版社，2003:1004-1007.。