二项Logistic回归参数最大似然估计的计算

二项Logistic回归参数最大似然估计的计算1Logistic回归的基本思想在地震风险评估中，研究者往往关心地震发生时，地表发生破裂的概率，地表破裂是由哪些因素引起的等问题。利用以往的相关数据找出统计规律性来解决这些问题，实质上可以转化为一个多元回归分析问题()Yfx，其中12[,,,]Tkxxxx，为随机变量。由于因变量Y的取值只有两个状态：破裂(1Y)和不破裂(0Y)，因此直接寻找因变量Y和自变量x的关系非常困难。于是，可以把研究问题转换一个角度，不去直接分析Y和x的关系，而是分析条件概率{1}PYx和x的关系，这等价于寻找一个取值在0到1之间的连续函数(){1}pPYxx。数学上满足这种条件的函数存在且不唯一，Logistic回归就是满足这种要求的函数之一。和线性回归分析类似，Logistic回归基本原理就是利用一组观测数据拟合一个Logistic模型，然后借助这个模型来揭示总体中若干个自变量与一个因变量取每个值的概率之间的依存关系，并评估用这一模型模拟相关事物变化规律的准确性。具体地说，Logistic回归分析可以从统计意义上确定在消除了其它变量的影响后，每一个自变量的变化是否引起因变量取某个值的概率的变化，并估计出在其它自变量固定不变的情况下，每个自变量对因变量取某个值的概率的数值影响大小。在使用Logistic回归分析前，需要明确模型的使用条件：1、要求因变量是分类变量，包括顺序变量和名义变量，不管哪种变量，都要用数字表示它，如可以令1Y表示地震发生时地表破裂，令0Y表示地震发生时地表未破裂；2、自变量可以是(i)数值型连续变量，如地震的震级，(ii)顺序变量，如覆盖层的厚度分组(10-20m，20-40m等)，(iii)名义变量，如地震类型，可令走滑型地震为1，正断型地震为2，逆冲型地震为3。2多元二项Logistic回归模型的定义由于地震发生时地表是否破裂受到多个因素的影响，故引入多元Logistic回归模型。假设因变量Y是一个取值为1和0的二值变量，12[,,,]Tkxxxx是影响Y的k个因素，回归系数01[,,,]Tkβ，则Y关于x的k元Logistic回归模型定义为0112201122exp()exp([1,])(){1}1exp()1exp([1,])TkkTkkxxxpPYxxxxβxxxβ(1)由式(1)可得1{0}1exp([1,])TPYxxβ(2)3Logistic回归参数估计我们用最大似然估计方法去求模型的参数。再假设从总体(,)Yx中抽取一个容量为12nn的随机样本112(1,),(1,),,(1,)nxxx，111212(0,),(0,),,(0,)nnnnxxx，其中12(,,,)iiiikxxxx，121,2,,inn，则有似然函数为112111211111011221101122011221(){1}{0}exp()11exp()1exp()exp([1,])11exp([1,])1exp([1,])nnniiiinnnniikikiiniikikiikiknTiTTiiniiLPYPYxxxxxxxxxβxxxβxβxβ121nn(3)两边取对数，整理可得112011221011221ln()()ln[1exp()]niikikinniikikiLxxxxxxβ(4)写成向量形式为11211ln()[1,]ln[1exp([1,])]nnnTTiiiLβxβxβ(4’)为求式(4)的驻点，可求对数似然函数ln()Lβ关于β的似然方程组为121211210111111001101110111111111011ln()exp()101exp()1exp()ln()exp()1exp()1nnnnikikiiikikikiknnnniikikiiiiiiikikLxxnnxxxxLxxxxxxxxββ12112112101101111110110110exp()ln()exp()01exp()1exp()nniikiknnnnnnikikikikikikiiiikikikikikxxLxxxxxxxxxxβ(5)写成向量形式为121121121101111111ln()101exp()ln()01exp()ln()01exp()nniinnniiiiinnnikikiikiLnLxxLxxβxββxββxβ式(5)为非线性方程组，一般情况下没有解析解，可以用Newton-Raphson迭代方法求其数值解，令12121121112110111111111011110111exp()11exp()1exp()()1ex1exp()nnknnijijjiinnnniiikiiijijijnnnikikkiijijjnxnxxxxxFxxxxββ121121110p()1exp()nniinnnikikiiixxxβxβ(6)则()Fβ关于β的Jacobian矩阵为1212120100111222111000111101exp()exp()exp()[1exp()][1exp()][1exp()]exp()()[1expkkkjijijijikjijnnnnnnjjjkkkiiijijjijjijjjjkijijjxxxxxxxxxxJβ12121212101011222111000111012101exp()exp()()][1exp()][1exp()]exp()[1exp()]kkijijiikjijnnnnnnjjkkkiiijijjijjijjjjkikjijnjkijijjxxxxxxxxxxx212121221001122110011121exp()exp()[1exp()][1exp()]exp()e[1exp()]kkikijijikjijnnnnnjjkkiijijjijjjnniiiixxxxxxxxxβxβ121212121222112111222111xp()exp()[1exp()][1exp()]exp()exp(exp()[1exp(][1exp()][1exp()]exp()[1ennnniikiiiiinnnnnniiiiiikiiiiiiiikixxxxxxxβxβxβxβxβxβxβxβxβxβxβ12121221222111exp()exp()xp()][1exp()][1exp()]nnnnnnikiiikiiiiiiixxxxβxβxβxβxβ(7)向量形式为12211exp()()1[1exp()]nniiTiiiJxββxxxβ(7’)根究Newton-Raphson方法的原理，可得参数β迭代公式为1(1)()()()()(),0,1,2,.nnnnJFnββββ(8)算法如下：Step1:给定参数β的初值参数(0)β和误差容许精度，令0n；Step2：计算1(1)()()()()(),0,1,2,.nnnnJFnββββ；Step3:若1()()()()nnJFββ或()()nFβ，即满足容许的精度，则结束，否则更新参数()(1)nnββ，1nn，转至Step2.当给定地震发生时，地表破裂是否发生的数据时，根据上面的算法，可以求出参数的最大似然估计。我们按照上述算法用MATLAB编写了多元Logistic回归参数估计的程序，可以给出参数估计值，详见附录。附录1用Newton-Raphson方法求解参数，附录2用直接优化对数似然函数方法求解参数，附录3用MATLAB自带的广义回归模型估计参数。附录4是别人做的Logistic回归的例子，用的软件是SAS（一种经过美国FDA认证的昂贵的商业统计软件）得到的结果。附录5是SPSS操作的过程和得到的结果。附录1:MatlabFilesforLogisticRegression%假设我们有一个数据，45个观测值，四个变量，包括：%1.age（年龄，数值型）；%2.vision（视力状况，分类型，1表示好，0表示有问题）；%3.drive（驾车教育，分类型，1表示参加过驾车教育，0表示没有）和%4.一个分类型输出变量accident（去年是否出过事故，1表示出过事故，0表示没有）。%我们的目的就是要考察前三个变量与发生事故的关系。%第1至4列分别为accidentagevisiondrive;clear,clc,closealldata=[117111440014810155001751103501042110570002801020010381004501047110520005501168101181016800148111170017011172101350111910162100391104011055000680102510017000450104401067000550116110119101690012311119001721117410131011161016110];Y=data(:,1);X=data(:,3:4);beta0=[0.110;1.7137;-1.5000]+1*rand(3,1);%rand(4,1);%猜测的初始值%自带的fsolve算法求解没有问题，核心原理也是Newton-Raphson算法%options=optimset('Display','iter','TolFun',1e-8)beta=fsolve(@(be)LogisticF(be,Y,X),beta0)%,options)%自编的Newton-Raphson算法，对初值比较敏感beta=LogisticRegressNR(Y,X,beta0)可以看到，自编函数与MATLAB自带函数Solve得到的结果相同，自编函数的缺点是对初值敏感，没有编写对应的策略，可用GA、PSO等算法定初始值。自编函数中用到的子函数：%%%子函数极大似然方程组functionF=LogisticF(beta0,Y,X)n=length(Y);%样本个数n=n1+n2X1=[ones(n,1)X];dims=size(X1,2);%待估参数个数ind1=(Y==1);%Y=1的样本个数col1=sum(X1(ind1,:))';%似然方程组F的第一部分col2=zeros(dims,1);%似然方程组F的第二部分初值%以下的向量表达好像不