二项Logistic回归参数最大似然估计的计算资料.docx
《二项Logistic回归参数最大似然估计的计算资料.docx》由会员分享，可在线阅读，更多相关《二项Logistic回归参数最大似然估计的计算资料.docx（27页珍藏版）》请在冰点文库上搜索。
二项Logistic回归参数最大似然估计的计算资料
二项Logistic回归参数最大似然估计的计算
1Logistic回归的基本思想
在地丧风险评估中，研究者往往关心地震发生时，地表发生破裂的概率，地表破裂是由哪些因素引起的等问题。
利用以往的相关数据找岀统讣规律性来解决这些问题，实质上可以转化为一个多元回归分析问题y=，其中召,农，…，兀r，£为随机变量。
由于因变量丫的取值只有两个状态：
破裂（y=i）和不破裂（r=o）,因此直接寻找因变量丫和自变量x的关系非常困难。
于是，可以把研究问题转换一个角度，不去直接分析Y和*的关系，而是分析条件概率p（r=i|x}和x的关系，这等价于寻找一个取值在o到1之间的连续函数P（x）=P{Y=l|x}o
数学上满足这种条件的函数存在且不唯一,Logistic回归就是满足这种要求的函数之一。
和线性回归分析类似.Logistic回归基本原理就是利用一组观测数据拟合一个Logistic模型,然后借助这个模型来揭示总体中若干个自变量与一个因变疑取每个值的概率之间的依存关系，并评估用这一模型模拟相关事物变化规律的准确性。
具体地说，Logistic回归分析可以从统计意义上确左在消除了其它变量的影响后，每一个自变量的变化是否引起因变量取某个值的概率的变化，并估计出在其它自变量固定不变的情况下，每个自变量对因变量取某个值的概率的数值影响大小。
在使用Logistic回归分析前，需要明确模型的使用条件：
1、要求因变量是分类变量，包括顺序变量和划义变量，不管哪种变量，都要用数字表示它，如可以令丫=1表示地震发生时地表破裂，令}7=0表示地丧发生时地表未破裂：
2、自变量可以是⑴数值型连续变量,如地丧的震级，（ii）顺序变量,如覆盖层的厚度分组（10~20m,20-40m等），（iii）名义变量,如地箴类型，可令走滑型地震为1,正断型地震为2,逆冲型地震为3。
2多元二项Logistic回归模型的定义
由于地震发生时地表是否破裂受到多个因素的影响，故引入多元Logistic回归模型。
假设因变量F是一个取值为1和0的二值变量，兀=[為宀，…,xj是影响丫的R个因素，回归系数“=[0。
胡,…,0J,则Y关于X的k元Logistic回归模型定义为心）=P{Y=i\x}=cxp（0（严阳+0壬+.・・+处丿=exp（[L八卩）
（1）
1+exp（仇+0內+/32x2+…+0內）1+exp（[1,x1\fl）
由式
（1）可得
P{Y=O\x}=T
1l+exp（［l,x7］fi）
3Logistic回归参数估计
我们用最大似然估计方法去求模型的参数。
再假设从总体（Y,x）中抽取一个容量为耳+”2的随机样本（1,旺）,（1山2），・・・，（1，乙”
（0心严1）,（0,%严2），・・・，（°心严），其中£=（兀］，兀・2,…心），j=12・・M+s则有似
然函数为
/1|眄+丿|,
£（0）=口川31|£}盯P（eo|£}
/-I
i■叫+1
=亓exp（炕+0』厂0N厂…+久$）1亓1⑶
V1+cxp（0o+朋|+02看2+…+0皿）去占1+exp（几+0丙+戸2齐2+…+以心）=计cxp（[l.打M），幵]
一*1+cxp（[Lx,r]0）七£1+exp（[Lxj|0）
两边取对数，整理可得
InUP）=力（0（）+0E\+际2+・•・+pkxik）
…⑷
一工ln[l+exp（A,+^xn+p2xi2+•••+pkxik）]
写成向量形式为
/tj+ru
讥（0）=£[1,x；]p-£ln[l+exp（[l,xT]fi）]（4J
为求式（4）的驻点，可求对数似然函数InL（fl）关于“的似然方程组为
=叫-》
/-I
=0
1+exp（-0（）-0®0"）
exp（0D+0]©+・・・+A.q）
1+exp（久+0心+…+0丹）
clnLip）
«•k\十也
=Z*Vh-Z
i・1r-1
心exp（0"+Q.D+・・・+0H』
1+exp（0。
+Q兀]+・・・+卩皿）
1-!
1-1
1+exp（-几一禹心
=0
0心
clnLip）
n,n,+n2
=Z-vrt-Si・lr-1
.Jexp（炕-0K：
+…+0"」_土,
l+exp（0|）+0].£]+・・・+0血）i.ini-i1+exp（-0仆一0[.勺0內）
写成向量形式为
卄乞—*—=0
800til+exp（-£0）
泌邑£补£儿=。
切tr台l+exp（-xQ
空込邑芝“=0
氓£山台l+exp（—*）
式（5）为非线性方程组，一般情况下没有解析解，可以用Ncwton-Raphson迭代方法求其数值解，令
F（"）=
Il+CXp（-0°-》0內）>-1
«!
坷.小丫
J—山I1+CXP（-仇一工0內）
MEI
nA_Y
匕l+cxpi”）
叫叫"：
r
yv-y亠：
zr11幺i+cxp（-x〃）
=0
m円+心r
x-^-z―—
I-l+CXp（-0°_2>A）
j"
则F（“）关于“的Jacobian矩阵为
5«1+心v
yv-y2ii
厶Atrl+expc-x,#）
i-i
J®=
…exp（-几-£几耳）
-Z——
z[l+exp（-0o-》M.%）]・
內杯•环exp（-几一》0声丿
；'
z[l+exp（—几一士0几）F
（
屮“、心exp（-0“一》Q.%）
-2r——
I[1+exp（-几一儿）r
4-1
1
屮„、時exp（-A,-內）z11+exp（-久-乞0儿）F
J-!
c
计.总exp（-几一》0心）
V」■】
•••/■
I[l+exp（m-》0r^）Fj-i
1
...
◎[1+exp（-鶴一乞0再厅j-i
入exp（-0。
-£0%）/-|
-Zx——
“[i+exp（-A.-ZM；）l2
A
n.林乙％exp（-炕一》0人）
“[1+exp（-久-乞0儿）F
-
j-i
导exp（-“0
它心exp（-“）台[l+exp（-x/）I2
它心exp（-x0Z?
[1+exp（-x,0）F
m|l+exp（-x/）]2
©exp（-兀0）纟[1+exp（-x（0F
£xJexp（-M0
Zf[l+exp（-x/）I2
_VAnA；Aexp（-x/>
Z>[l+exp（-x.0）F
W.和exp（-x,0）
—旨从•环exp（-x,“）
它時exp（-〃）台[l+exp（-x/）F
[l+exp（-x,/?
）r
m[l+expC-x/）]2
屮与•叮°"<-几一乞0.4）
亍2
•••—7，.
I[l+eKp（m£Wj-i
（7）
向量形式为
根究Newton-Raphson方法的原理，可得参数“迭代公式为
卩e=0"-[丿0”F（/T））丿=0丄2,….⑻
算法如下：
Step1:
给左参数“的初值参数“®和误差容许精度s令〃=0：
sg2：
计算“曲）=/r>—p（/r））『F（/r）M=oj2…•：
Step3:
若
V£，即满足容许的精度，则结束，否
则更新参数p^=p{n^\n=n+\,转至Stcp2・
当给左地後发生时，地表破裂是否发生的数据时，根据上而的算法，可以求出参数的最大似然估计。
我们按照上述算法用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.一个分类型输出变Maccident（去年是否岀过事故，1表示岀过事故，0表示没
有）。
%我们的目的就是要考察前三个变量与发生事故的关系。
%第1至4列分别为accidentagevisiondrive;
cleanclc,closeall
data=[11711
14400
14810
15500
17511
03501
04211
05700
02801
02001
038100450104711
05200
05501
16810
11810
16800
14811
11700
17011
17210
13501
11910
16210
03911
04011
05500
06801
02510
01700
04501
04401
06700
05501
16110
11910
16900
12311
11900
17211
17410
13101
11610
16110];
Y=data（:
J）;
X=data（:
3:
4）;
beta0=[0.110;1・7137;-1.5000]+1*rand（3,1）;%rand（4,1）;%猜测的初始值
%自带的fsolvc算法求解没有问题，核心原理也是Newton-Raphson算法
%options=optimset（,Display\'iter\TolFun\1e-8）
beta=fsoIve（@（be）LogisticF（be,Y,X）,betaO）%,options）
%自编的Newton-Raphson算法，对初值比较敏感
beta=LogisticRcgressNR（Y.X.betaO）
可以看到，自编函数与MATLAB自带函数Solve得到的结果相同，自编函数的缺点是对初值敏感，没有编写对应的策略，可用GA、PSO等算法泄初始值。
自编函数中用到的子函数：
%%%子函数极大似然方程组
functionF=LogisticF（betaO,Y,X）
n=length（Y）;%样本个数n=nl+n2
XI=[ones（nj）X];
dims=size（Xl,2）;%待估参数个数
indl=（Y==l）;%Y=1的样本个数
coll=sum（XKindi,:
））';%似然方程组F的第一部分
col2=zeros（dims,l）;%似然方程组F的第二部分初值
%以下的向量表达好像不对
%col2=sum（X1./（1+exp（X1♦betaO）））1;
%
fori=l:
n
col2=col2+（Xl（i,:
）/（l+exp（-Xl（i,:
）*betaO））y；
end
F=[coll-co!
2J;
%Newton-Rapson算法中的Jacobian矩阵
functionM=LogisticJM（betaO,Y,X）
n=length（Y）;%样本个数
XI=[ones（nJ）X];%变量个数
dims=size（Xl,2）;
M=zeros（dims）;
fori=l:
n
M=M+Xl（i,:
）,*Xl（i,:
）*exp（-X1（i,:
）*betaO）/（（1+exp（-X1（i,:
）*betaO））A2）;end
M=-M;
functionbeta=LogisticRcgressNR（YX.bctaO）
%用牛顿一拉普生方法求极大似然估计
%Y因变量样本观测值，取值为1,表示事件发生，取值为0,表示事件不发生
%X多个自变量的样本观测值，默认X的第一列全为1
%betaO猜测的beta的初始值
%%%%王福昌2015-6-15
%%主函数
itermax=10;%最多迭代次数
errstol=le-4;%容忍误差
iters=0;%迭代次数
%[n.k]=size（X）;诙n样本容量，k自变量个数deltabeta=ones（size（betaO））;
betal=betaO+deltabeta;%虚拟迭代误差向量
while（iterserrstol）
deltabeta=・LogisticJM（betaO,Y.X）\LogisticF（betaO,Y,X）;
betal=betaO+deltabeta;
betaO=beta1;iters=iters+1;
end
beta=betaO;
附录2:
直接优化对数似然函数，也能得到同样结果
functionF=LogisticRegressOpt（betarYX）
%用最优化方法求极大似然估计
%Y因变量样本观测值，取值为1,表示事件发生，取值为0,表示事件不发生
%X多个自变量的样本观测值，默认X的第一列全为1
%betaO猜测的beta的初始值
%%%%王福昌2015-6-15
%%主函数
%迭代次数
%极大似然函数
indl=（Y==l）;
n=length（Y）;
XI=[ones（nj）X];
F=sum（X1（ind1,:
）*beta）-sum（Iog（l+exp（Xl*beta）））;
F=-F;
%假设我们有一个数据，45个观测值，四个变量，包括:
%1.
%2.
%3.
%4.有）。
age（年龄,数值型）；
vision（视力状况，分类型，1表示好，0表示有问题）：
drive（驾车教冇，分类型，1表示参加过驾车教育，0表示没有）和
一个分类型输出变Maccident（去年是否出过事故，1表示出过事故，0表示没
%我们的目的就是要考察前三个变量与发生事故的关系。
%第1至4列分别为accidentagevisiondrive;
data=[11711
14400
数据与附录1中相同
11610
16110];
Y=data（:
l）;
X=data（:
3:
4）;
beta0=[l;l;0];%[0.1110;1.7137;-1.5000];%rand（4J）;%猎测的初始值
Ibetafval]=fminsearch（@（beta）LogisticRcgressOpt（beta,Y.X）,bctaO）
可见，参数估计结果相同
附录3:
b=glmfit（X,Y,'binomial*,link*,logit*）
p=glmval（beta,X,logit'）
可以看到，与前面的两种方法结果相同。
附录4：
对比的例子
proclogistic（logistic回归的SAS实现一无哑变量）
（2012-03-1321:
30:
54）
转拔▼
原文mt：
proclogistic（logisticH归的SAS实现一无呵变放）作者，贝搭数据统i|•工作
Logistic阿归主耍用来处別应变肚为分类变煲的何題•比如生存和死亡•患病和未忠擒乐当然研宪者咲心的问题是哪牛因索导致了患釣或不忠衲.哪些足生存和死亡.
Logistic的sas语句很简单.其基本语句见下：
PROCLOGISTICDATA=SAS-data-set;
MODELresponse=independents;
BYvariables;
OUTPUT
>:
WEIGHTvariable:
Proclogistic语句戏认计味应变址值眾小｛阴性结果一股赋值为0）的概率.但足常密我们想要得到的是阳性給果的概率•即
赋的数值的概率（二分类变fit时一般赋値为1）.选项Mdescending解决了这一何题。
Model语句用干定义应变贵和自变址：
实例：
锻设我们有一个数据•45个观测值.四个变fit.包括:
age（年龄.数值型〉：
2.vision（视力状况.分类型.1表示好.0茨示有何题〉：
3.drive（驾车教育，分类空.1表示参加过驾车教育•0表示没有）和
4.一个分类製输出变®accident（去年足否出过审故.1表示岀过爭故.0茨示没有）。
我们的目的就足要考察前：
个变
壇与发生事故的关系。
datalogistic;
inputaccidentagevisiondrive:
datalines:
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
04501
04401
06700
05501
16110
11910
16900
12311
11900
17211
17410
13101
11610
16110run;
（1）proclogisticdata=logisticdescending;
modelaccident=agevisiondrive：
run;
如果想要在选择适当的自变fit筛选方法則使用一下语句：
（2）proclogisticdata=logisticdescending;
modelaccident=agevisiondrive/selection=：
stepwisesle=0.15sls=0.15stb:
run:
Selection用于选择筛选自变就的方法.有backward（向后法）、forward（向前法）、stepwise（逐步法）、score（最优子集法）、none（完全法〉五个选项.戏认为none：
SLE=概率值.入选标准.規定变掀入选怏型的显著性水平.前进法的炊认足0.5・逐步浓足0.15
SLS=ftt率值.剔除标准.指定变扯保留在怏里的显著水平.后退法炊认为0.10.逐步浓是0.15标准化偏回归系数STB可用来比较备个自变fit作用的人小
还可以输出贾信区何.语句如“
（3）proclogisticdata=logisticdescending;
modelaccident=agevisiondrive/selection=stepwisesle=0.15sls=0.15stbcl;run;
结呆：
TheLOGISTICProcedure
ModelInformation
DataSetWORK.LOGISTIC
ResponseVariable
acddent
NumberofResponseLevels
2
NumberofObservations
45
Model
binarylogit
OptimizationTechnique
Fishersscoring
ResponseProfile
OrderedValueaccident
TotalFrequency
1125
2020
Probabilitymodeledisaccident=d・
⑴给岀了木枳电的慕本信息•总思人釦‘I明。
缶要注您的足ResponseProfile中.accident=1排在忤位。
帕而我们说过.SAS的Logistic回归方楼log（odds）IV：
认的形式是处典那个变磧（ft比较小的・加J:
descending选项后・accident=1就挣在闫位了.
（2）ForwardSelectionProcedure
Step0.Interceptentered:
ModelConvergenceStatus
Convergencecriterion（GCONV=1E-8）satisfied.
ResidualChi-SquareTest
Chi-SquareDF
Pr>ChiSq
10.70573
0.0134
Step1・Effectvisionentered:
ModelConvergenceStatus
Convergencecriterion（GCONV=1E-8）satisfied.
ModelFitStatistics
CriterionInterceptOnlyInterceptandCovariates
SC65.63362.857
-2LogL61.82755.244
TestingGlobalNullHypothesis:
BETA=0
Test
Chi-Square
DF
Pr>ChiSq
LikelihoodRatio
6.5830
1
0.0103
Score
6.4209
1
0.0113
Wald
6.0756
1
0.0137
ResidualChi-SquareTest
Chi-Square
DFPr>ChiSq
4.981
升级会员 