二项Logistic 回归参数最大似然估计的计算资料.docx
- 文档编号:2050138
- 上传时间:2022-10-26
- 格式:DOCX
- 页数:16
- 大小:481.51KB
二项Logistic 回归参数最大似然估计的计算资料.docx
《二项Logistic 回归参数最大似然估计的计算资料.docx》由会员分享,可在线阅读,更多相关《二项Logistic 回归参数最大似然估计的计算资料.docx(16页珍藏版)》请在冰豆网上搜索。
二项Logistic回归参数最大似然估计的计算资料
二项Logistic回归参数最大似然估计的计算
1Logistic回归的基本思想
在地震风险评估中,研究者往往关心地震发生时,地表发生破裂的概率,地表破裂是由哪些因素引起的等问题。
利用以往的相关数据找出统计规律性来解决这些问题,实质上可以转化为一个多元回归分析问题,其中,为随机变量。
由于因变量的取值只有两个状态:
破裂()和不破裂(),因此直接寻找因变量和自变量的关系非常困难。
于是,可以把研究问题转换一个角度,不去直接分析和的关系,而是分析条件概率和的关系,这等价于寻找一个取值在0到1之间的连续函数。
数学上满足这种条件的函数存在且不唯一,Logistic回归就是满足这种要求的函数之一。
和线性回归分析类似,Logistic回归基本原理就是利用一组观测数据拟合一个Logistic模型,然后借助这个模型来揭示总体中若干个自变量与一个因变量取每个值的概率之间的依存关系,并评估用这一模型模拟相关事物变化规律的准确性。
具体地说,Logistic回归分析可以从统计意义上确定在消除了其它变量的影响后,每一个自变量的变化是否引起因变量取某个值的概率的变化,并估计出在其它自变量固定不变的情况下,每个自变量对因变量取某个值的概率的数值影响大小。
在使用Logistic回归分析前,需要明确模型的使用条件:
1、要求因变量是分类变量,包括顺序变量和名义变量,不管哪种变量,都要用数字表示它,如可以令表示地震发生时地表破裂,令表示地震发生时地表未破裂;2、自变量可以是(i)数值型连续变量,如地震的震级,(ii)顺序变量,如覆盖层的厚度分组(10-20m,20-40m等),(iii)名义变量,如地震类型,可令走滑型地震为1,正断型地震为2,逆冲型地震为3。
2多元二项Logistic回归模型的定义
由于地震发生时地表是否破裂受到多个因素的影响,故引入多元Logistic回归模型。
假设因变量是一个取值为1和0的二值变量,是影响的个因素,回归系数,则关于的元Logistic回归模型定义为
(1)
由式
(1)可得
(2)
3Logistic回归参数估计
我们用最大似然估计方法去求模型的参数。
再假设从总体中抽取一个容量为的随机样本,,其中,,则有似然函数为
(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
04501
04401
06700
05501
16110
11910
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)
n=length(Y);%样本个数
X1=[ones(n,1)X];%变量个数
dims=size(X1,2);
M=zeros(dims);
fori=1:
n
M=M+X1(i,:
)'*X1(i,:
)*exp(-X1(i,:
)*beta0)/((1+exp(-X1(i,:
)*beta0))^2);
end
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
deltabeta=-LogisticJM(beta0,Y,X)\LogisticF(beta0,Y,X);
beta1=beta0+deltabeta;
beta0=beta1;iters=iters+1;
end
beta=beta0;
附录2:
直接优化对数似然函数,也能得到同样结果
functionF=LogisticRegressOpt(beta,Y,X)
%用最优化方法求极大似然估计
%Y因变量样本观测值,取值为1,表示事件发生,取值为0,表示事件不发生
%X多个自变量的样本观测值,默认X的第一列全为1
%beta0猜测的beta的初始值
%%%%王福昌2015-6-15
%%主函数
%迭代次数
%极大似然函数
ind1=(Y==1);
n=length(Y);
X1=[ones(n,1)X];
F=sum(X1(ind1,:
)*beta)-sum(log(1+exp(X1*beta)));
F=-F;
%假设我们有一个数据,45个观测值,四个变量,包括:
%1.age(年龄,数值型);
%2.vision(视力状况,分类型,1表示好,0表示有问题);
%3.drive(驾车教育,分类型,1表示参加过驾车教育,0表示没有)和
%4.一个分类型输出变量accident(去年是否出过事故,1表示出过事故,0表示没有)。
%我们的目的就是要考察前三个变量与发生事故的关系。
%第1至4列分别为accidentagevisiondrive;
data=[11711
14400
…………数据与附录1中相同
11610
16110];
Y=data(:
1);
X=data(:
3:
4);
beta0=[1;1;0];%[0.1110;1.7137;-1.5000];%rand(4,1);%猜测的初始值
[betafval]=fminsearch(@(beta)LogisticRegressOpt(beta,Y,X),beta0)
可见,参数估计结果相同
附录3:
b=glmfit(X,Y,'binomial','link','logit')
p=glmval(beta,X,'logit')
可以看到,与前面的两种方法结果相同。
附录4:
对比的例子
proc logistic(logistic回归的SAS实现--无哑变量)
(2012-03-1321:
30:
54)
转载▼
原文地址:
proc logistic(logistic回归的SAS实现--无哑变量
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 二项Logistic 回归参数最大似然估计的计算资料 Logistic 回归 参数 最大 估计 计算 资料