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

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

1 Possion回归模型的定义
假设因变量是一个服从Poission分布的随机变量，是影响的个因素，是协变量向量，是回归参数向量，则关于的元Poission回归模型定义为
(1)
其中
2 参数估计
我们用最大似然估计方法去求模型的参数。
假设从总体中抽取一个容量为的随机样本，其中，则有似然函数为
(2)
两边取对数，整理可得
(3)
为研究方便，以下不妨记。为求式(3)的最大值点，即最大似然估计，可求对数似然函数关于的似然方程组为
， (4)
具体形式为
(5)
式(5) 为非线性方程组，一般情况下没有解析解，可以用Newton-Raphson迭代方法求其数值解，令
(6)
则关于的Jacobian矩阵为
(7)
具体形式为
(7)
对应的向量形式为
(7’)
根据Newton-Raphson方法的原理，可得参数迭代公式为
(8)
算法如下：
Step 1: 给定参数的初值参数和误差容许精度，令；
Step 2：计算；
Step 3: 若，即满足容许的精度，则结束，否则更新参数，，转至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.