二项Logistic 回归参数最大似然估计的计算资料Word格式.docx

1、再假设从总体中抽取一个容量为的随机样本，其中，则有似然函数为 （3）两边取对数，整理可得 （4）写成向量形式为 （4）为求式（4）的驻点，可求对数似然函数关于的似然方程组为 （5）式（5） 为非线性方程组，一般情况下没有解析解，可以用Newton-Raphson迭代方法求其数值解，令 （6）则关于的Jacobian矩阵为 （7）向量形式为 （7）根究Newton-Raphson方法的原理，可得参数迭代公式为 （8）算法如下：Step 1: 给定参数的初值参数和误差容许精度，令；Step 2：计算；Step 3: 若或，即满足容许的精度，则结束，否则更新参数，转至Step2.当给定地震发生时，地

2、表破裂是否发生的数据时，根据上面的算法，可以求出参数的最大似然估计。我们按照上述算法用MATLAB编写了多元Logistic回归参数估计的程序，可以给出参数估计值，详见附录。附录1用Newton-Raphson方法求解参数，附录2用直接优化对数似然函数方法求解参数，附录3用MATLAB自带的广义回归模型估计参数。附录4是别人做的Logistic回归的例子，用的软件是SAS（一种经过美国FDA认证的昂贵的商业统计软件）得到的结果。附录5是SPSS操作的过程和得到的结果。附录1:Matlab Files for Logistic Regression% 假设我们有一个数据，45个观测值，四个变量，

3、包括：% 1. age（年龄，数值型）；% 2. vision（视力状况，分类型，1表示好，0表示有问题）；% 3. drive（驾车教育，分类型，1表示参加过驾车教育，0表示没有） 和% 4. 一个分类型输出变量accident（去年是否出过事故，1表示出过事故，0表示没有）。% 我们的目的就是要考察前三个变量与发生事故的关系。 % 第1至4列分别为 accident age vision drive ; clear,clc,close alldata = 1 17 1 11 44 0 01 48 1 01 55 0 01 75 1 10 35 0 10 42 1 10 57 0 00 28

4、 0 10 20 0 10 38 1 00 45 0 10 47 1 10 52 0 00 55 0 11 68 1 01 18 1 01 68 0 01 48 1 11 17 0 01 70 1 11 72 1 01 35 0 11 19 1 01 62 1 00 39 1 10 40 1 10 55 0 00 68 0 10 25 1 00 17 0 00 44 0 10 67 0 01 61 1 01 69 0 01 23 1 11 19 0 01 72 1 11 74 1 01 31 0 11 16 1 01 61 1 0;Y = data（:,1）;X= data（:,3:4）;be

5、ta0 = 0.110 ;1.7137 ; -1.5000+1*rand（3,1）;%rand（4,1）; %猜测的初始值%自带的fsolve 算法求解没有问题，核心原理也是Newton-Raphson 算法%options = optimset（Display,iterTolFun,1e-8）beta = fsolve（be）LogisticF（be,Y,X）,beta0）%,options） % 自编的Newton-Raphson 算法，对初值比较敏感beta = LogisticRegressNR（Y,X,beta0）可以看到，自编函数与MATLAB自带函数Solve得到的结果相同，自编

6、函数的缺点是对初值敏感，没有编写对应的策略，可用GA、PSO等算法定初始值。自编函数中用到的子函数：%子函数 极大似然方程组function F = 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 的第二部分初值% 以下的向量表达好像不对% col2 = sum（X1

7、./（1+exp（X1*beta0）% for i = 1:n col2 = col2 + （X1（i,:）/（1+exp（-X1（i,:）*beta0）endF = col1 - col2;% Newton-Rapson算法中的 Jacobian 矩阵function M = 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;function beta = LogisticRegre

8、ssNR（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;% 虚拟迭代误

9、差向量while （iterserrstol） deltabeta = -LogisticJM（beta0,Y,X）LogisticF（beta0,Y,X）; beta1 = beta0 + deltabeta; beta0 = beta1; iters = iters +1;beta = beta0;附录2: 直接优化对数似然函数，也能得到同样结果function F = LogisticRegressOpt（beta,Y,X）% 用最优化方法求极大似然估计% 迭代次数% 极大似然函数ind1 = （Y=1）;F = sum（X1（ind1,:）*beta） - sum（log（ 1+exp

10、（X1*beta）;F = -F;数据与附录1中相同beta0 = 1;1;0;% 0.1110 ; -1.5000;beta fval = fminsearch（beta） LogisticRegressOpt（beta,Y,X）,beta0）可见，参数估计结果相同附录3:b =glmfit（X,Y,binomial, linklogit）p = glmval（beta,X, 可以看到，与前面的两种方法结果相同。附录4：对比的例子proclogistic（logistic回归的SAS实现-无哑变量）（2012-03-13 21:30:54）转载原文地址：logistic（logistic回归的SAS实现-无哑变量