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

Poisson回归参数最大似然估计的计算1Possion回归模型的定义假设因变量Y是一个服从Poission分布的随机变量，12,,,qxxx是影响Y的k个因素，12[1,,,,]TqXxxx是协变量向量，01[,,,]Tqβ是回归参数向量，则Y关于x的k元Poission回归模型定义为(){}exp(()),0,1,2,.!kXPYkXXkk(1)其中()exp()0.TXXβ2参数估计我们用最大似然估计方法去求模型的参数。假设从总体(,)YX中抽取一个容量为n的随机样本1122(,),(,),,(,)nnyXyXyX，其中12[1,,,,],1,2,,TkkkkqXxxxkn，则有似然函数为11exp()(){}exp(exp())!kyTnnTkkkkkkkXLPYyXXyβββ(2)两边取对数，整理可得1ln()exp()ln(!)nTTkkkkkLyXXyβββ(3)为研究方便，以下不妨记01kx。为求式(3)的最大值点，即最大似然估计，可求对数似然函数ln()Lβ关于β的似然方程组为1ln()[exp()]nTkkikikkiLyxxXββ，0,1,,.iq(4)具体形式为1011111ln()[exp()]0ln()[exp()]0ln()[exp()]0nTkkknTkkkkknTkkqkqkkkLyXLyxxXLyxxXββββββ(5)式(5)为非线性方程组，一般情况下没有解析解，可以用Newton-Raphson迭代方法求其数值解，令11111[exp()][exp()]()[exp()]nTkkknTkkkkknTkkqkqkkyXyxxXFyxxXββββ(6)则()Fβ关于β的Jacobian矩阵为21ln()()exp(),0,1,,,0,1,,.nTkikjkijLJxxXiqjqββ(7)具体形式为1111211111121111exp()exp()exp()exp()exp()exp()()exp()exp()exp()nnnTTTkkkkqkkkknnnTTTkkkkkkqkkkknnnTTTkqkkqkkkqkkkkXxXxXxXxXxxXJxXxxXxXββββββββββ(7)对应的向量形式为1()exp()nTTkkkkJXXXββ(7’)根据Newton-Raphson方法的原理，可得参数β迭代公式为1(1)()()()()(),0,1,2,.mmmmJFmββββ(8)算法如下：Step1:给定参数β的初值参数(0)β和误差容许精度，令0m；Step2：计算1(1)()()()()(),0,1,2,.mmmmJFmββββ；Step3:若()()mFβ，即满足容许的精度，则结束，否则更新参数()(1)mmββ，1mm，转至Step2.functionF=PoissionRegressopt(b,Y,X)n=length(Y);F=0;fork=1:nF=F+Y(k)*X(k,:)*b-exp(X(k,:)*b);%-factorial(Y(k));endF=-F;functionF=PoissionF(b,Y,X)n=length(Y);F=zeros(size(b));fork=1:nF=F+Y(k)*X(k,:)'-exp(X(k,:)*b)*X(k,:)';endfunctionJM=PoissionJM(b,Y,X)n=length(Y);JM=zeros(size(b,1));fork=1:nJM=JM+exp(X(k,:)*b)*X(k,:)'*X(k,:);endfunction[bmfv1,fv2]=PoissionNR(bm0,Y,X)itermax=30;errstol=1e-4;iters=0;deltabm=ones(size(bm0));bm1=bm0+deltabm;while(itersitermax)||(max(abs(deltabm))errstol)deltabm=pinv(PoissionJM(bm0,Y,X))*PoissionF(bm0,Y,X);bm1=bm0+deltabm;bm0=bm1;iters=iters+1;endbm=bm0;fv1=PoissionF(bm,Y,X);fv2=PoissionRegressopt(bm,Y,X);附录1:b=glmfit(X0,Y,'poisson','log')b=1.50430.45180.35780.2388可以看到，结果一致。比文献【1】中的结果要好一点参考文献【1】茆诗松主编.统计手册[M].北京:科学出版社，2003:1004-1007.