二项Logistic 回归参数最大似然估计的计算资料Word格式.docx
《二项Logistic 回归参数最大似然估计的计算资料Word格式.docx》由会员分享，可在线阅读，更多相关《二项Logistic 回归参数最大似然估计的计算资料Word格式.docx（16页珍藏版）》请在冰豆网上搜索。
再假设从总体中抽取一个容量为的随机样本，，其中，，则有似然函数为
（3）
两边取对数，整理可得
（4）
写成向量形式为
（4’）
为求式（4）的驻点，可求对数似然函数关于的似然方程组为
（5）
式（5）为非线性方程组，一般情况下没有解析解，可以用Newton-Raphson迭代方法求其数值解，令
（6）
则关于的Jacobian矩阵为
（7）
向量形式为
（7’）
根究Newton-Raphson方法的原理，可得参数迭代公式为
（8）
算法如下：
Step1:
给定参数的初值参数和误差容许精度，令；
Step2：
计算；
Step3:
若或，即满足容许的精度，则结束，否则更新参数，，转至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,closeall
data=[11711
14400
14810
15500
17511
03501
04211
05700
02801
02001
03810
04501
04711
05200
05501
16810
11810
16800
14811
11700
17011
17210
13501
11910
16210
03911
04011
05500
06801
02510
01700
04401
06700
16110
16900
12311
11900
17211
17410
13101
11610
16110];
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+n2
X1=[ones（n,1）X];
dims=size（X1,2）;
%待估参数个数
ind1=（Y==1）;
%Y=1的样本个数
col1=sum（X1（ind1,:
））'
;
%似然方程组F的第一部分
col2=zeros（dims,1）;
%似然方程组F的第二部分初值
%以下的向量表达好像不对
%col2=sum（X1./（1+exp（X1*beta0）））'
%
fori=1:
n
col2=col2+（X1（i,:
）/（1+exp（-X1（i,:
）*beta0）））'
end
F=[col1-col2];
%Newton-Rapson算法中的Jacobian矩阵
functionM=LogisticJM（beta0,Y,X）
%样本个数
%变量个数
M=zeros（dims）;
M=M+X1（i,:
）'
*X1（i,:
）*exp（-X1（i,:
）*beta0）/（（1+exp（-X1（i,:
）*beta0））^2）;
M=-M;
functionbeta=LogisticRegressNR（Y,X,beta0）
%用牛顿—拉普生方法求极大似然估计
%Y因变量样本观测值，取值为1，表示事件发生，取值为0，表示事件不发生
%X多个自变量的样本观测值,默认X的第一列全为1
%beta0猜测的beta的初始值
%%%%王福昌2015-6-15
%%主函数
itermax=10;
%最多迭代次数
errstol=1e-4;
%容忍误差
iters=0;
%迭代次数
%[n,k]=size（X）;
%n样本容量，k自变量个数
deltabeta=ones（size（beta0））;
beta1=beta0+deltabeta;
%虚拟迭代误差向量
while（iters<
itermax）||（max（abs（deltabeta））>
errstol）
deltabeta=-LogisticJM（beta0,Y,X）\LogisticF（beta0,Y,X）;
beta1=beta0+deltabeta;
beta0=beta1;
iters=iters+1;
beta=beta0;
附录2:
直接优化对数似然函数，也能得到同样结果
functionF=LogisticRegressOpt（beta,Y,X）
%用最优化方法求极大似然估计
%迭代次数
%极大似然函数
ind1=（Y==1）;
F=sum（X1（ind1,:
）*beta）-sum（log（1+exp（X1*beta）））;
F=-F;
…………数据与附录1中相同
beta0=[1;
1;
0];
%[0.1110;
-1.5000];
[betafval]=fminsearch（@（beta）LogisticRegressOpt（beta,Y,X）,beta0）
可见，参数估计结果相同
附录3:
b=glmfit（X,Y,'
binomial'
'
link'
logit'
）
p=glmval（beta,X,'
可以看到，与前面的两种方法结果相同。
附录4：
对比的例子
proc
logistic（logistic回归的SAS实现--无哑变量）
（2012-03-1321:
30:
54）
转载▼
原文地址：
logistic（logistic回归的SAS实现--无哑变量
升级会员 