印明昂, 王钰烁, 孙志礼, 郭兵
东北大学 机械工程与自动化学院，辽宁 沈阳 110819
收稿日期：2018-12-26
基金项目：国家自然科学基金资助项目(51775097，51875095)。
作者简介：印明昂(1985-), 男, 辽宁沈阳人, 东北大学讲师, 博士;
孙志礼(1957-), 男, 山东巨野人, 东北大学教授, 博士生导师。

摘要：基于Newton迭代法提出一种累积Logistic回归模型的参数估计方法.分析了迭代初值选取、常系数大小关系以及迭代过程Hessian矩阵奇异等影响算法收敛性的主要问题.通过自适应地选取迭代初值和控制迭代过程，避免了Hessian矩阵奇异的情况.利用美国凯斯西储大学轴承数据库(CWRU)数据进行验证，实验结果表明，本文方法在模型训练和验证的准确率上均高于统计学软件SPSS.并利用Booststrap随机试验验证了所提算法的稳健性.
关键词：累积Logistic回归参数估计Hessian矩阵自适应算法滚动轴承健康状态评估
Parameter Estimation of the Cumulative Logistic Regression Model Based on the Newton Iteration Method
YIN Ming-ang, WANG Yu-shuo, SUN Zhi-li, GUO Bing
School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
Corresponding author: YIN Ming-ang, E-mail: yinma@mail.neu.edu.cn.

Abstract: Based on the Newton iteration method, a parameter estimation method of cumulative Logistic regression was proposed. The selection of primary values of iteration, the relationship of sizes of constant coefficients and the singularity of Hessian matrix were analyzed, which are the main problems influencing stypticity in the process of iteration. Through the self-adaptive selection of primary values of iteration and control of the iteration process, the singularity of Hessian matrix was avoided. The data of the bearing database of Case Western Reserve University was used for validation. The experiment results show that the accuracy of model training and verification proposed is higher than that of the statistics software SPSS. Moreover, the robustness of the proposed algorithm is proved by the Booststrap random testing method.
Key words: cumulative Logistic regressionparameter estimationHessian matrixadaptive algorithmrolling bearinghealth state assessment
累积Logistic回归又称为有序多分类Logistic回归、次序Logistic回归^[1]，其具有模型变量少、分类效率高的特点.近年来，该方法在医药、金融及制造等领域得到了广泛的应用.文献[2]选取2015年3月至2017年4月接受培非格司亭(Pegfilgrastim)联合化疗的166例癌症案例，采用多变量有序Logistic回归分析，确定符合Pegfilgrastim预防以维持相对剂量强度的预测因素.文献[3]基于发动机健康变化的Logistic回归模型，实现了贝叶斯状态估计的健康状态序列更新，进而预测引擎未来的健康变化.文献[4]通过审查2003至2011年外来资本流入的综合影响，检验了次序Logistic回归模型应用于基于柯布-道格拉斯生产函数的实证经济增长研究的可行性.
在应用领域日益扩展，数据数量和维度不断增大的现状下，目前，累积Logistic回归模型的参数估计通常采用统计学软件SPSS或SAS的自带功能实现.本文基于Newton迭代法提出一种累积Logistic回归模型的参数估计方法，对其中比较敏感的迭代初值选取及迭代过程Hessian矩阵奇异问题进行了自适应控制，最后利用美国凯斯西储大学轴承数据库(CWRU)^[5]数据对该算法模型与SPSS模型进行了对比验证.
1 参数估计1.1 迭代格式的构建设(x, y)为有序J分类样本，其中输入观测值x为n×M矩阵，n代表容量，M代表维度；有序输出观测值y为n×1列向量，其元素y_i, i∈{1, 2, …, n}, y_i可取值为1, 2, …, J.则第i个输入x_i=[x_i1, x_i2, …, x_iM]的前j(j∈{1, 2, …, J})类累积Logistic表达式为^[6]

(1)

其中：p(y_i≤j|x_i)为累积概率(简记为p_ij)，且有p(y_i≤J|x_i)=1；α_j为截距，对于累积Logistic模型共有J-1个；β_m为回归系数，所有分类的系数相同.由此可得

(2)

则该输入属于第j类的概率为

(3)

最后由概率最大者确定其所属类别.不妨设Y为n×J矩阵，其元素y_ij为

(4)

则它的对数似然函数可表示为

(5)

由Newton迭代公式得到迭代格式：

(6)

其中，b_k=[α_k^T, β_k^T]^T=[α_1k, …, α_(J-1)k, β_1k, …, β_Mk]^T为未知参数的第k次迭代值，k=0, 1, 2, …；H(b_k)，f(b_k)分别为

(7)

(8)

H(b_k)元素如下

(9)

(10)

(11)

其中，r, j=1, 2, …, J-1;s, m=1, 2, …, M.
1.2 迭代初值的选取Newton迭代法虽然具有平方收敛速度，但它对迭代初值非常敏感，如果选取不当往往会导致迭代发散.下面通过数学推导，给出一种初值的选取方法.
首先将数据归一化，使所有输入数据统一映射到[0, 1]区间上^[7], 即

(12)

其中：x_im表示输入i第m列的数据；x_m为所有第m列数据组成的向量，i=1, 2, …, n，m=1, 2, …, M；min(·)，max(·)分别表示向量中的最小、最大值.
根据Fisher判别思想，通常各个分类中心差距明显^[8].以各类中的样本均值表示分类中心，如j类中心x_j=[x_j1，x_j2，…，x_jM]为所有属于第j类输入的样本均值.计算得出各类中心X=[x₁^T, x₂^T, …, x_J^T]^T作为选取初值的特定数据.
设初值为b₀=[α₀^T, β₀^T]^T，α₀=[α₁₀, …, α_(J-1)0]^T为常系数初值，β₀=[β₁₀, …, β_M0]^T为变量系数初值.将X和b₀代入式(2)，式(3)，得到以未知参数表示的各分类条件预测概率矩阵Q：

(13)

其中，Q_ij=P(y=j|x_i, b₀).根据定义，预测条件概率Q_jj应为所在行的最大值，显然，Q中只有三主对角线上的元素显著不为零.由此可得下列J个关于b₀中各元素的不等式：

(14)

为使x_j属于第j类的概率最大，推荐取p₁=0.3，p₂=0.7.同时，为避免β₀各元素间差距过大，导致Hessian矩阵奇异，利用系数f_d1>0对其进行限定：

(15)

由式(14)可推出：

(16)

其中，c=ln(1/p₁-1) -ln(1/p₂-1).进而有

(17)

式(14)~式(17)给出了b₀的大致范围.为进一步压缩区间，提高选取成功率，将b₀分为三部分：b₀=[α₀^T, β₀₁^T, β₀₂^T]^T，其中α₀=[α₁₀, …, α_(J-1)0]^T，β₀₁=[β₁₀, …, β_k0]^T，β₀₂=[β_(k+1)0, …, β_M0]^T，k∈[1, J].将式(14)的第一和最后一式进行缩放：

(18)

(19)

其中f_d2=M·f_d1.设α₀，β₀₁已知，代入式(14)整理为关于β₀₂的不等式组：

(20)

其中，c₁=ln(1/p₁-1)，c₂=ln(1/p₂-1)；T_i1=[x_i1, x_i2, …, x_ik]，T_i2=[x_i(k+1), …, x_iM]，i=1, 2, …, J.将式(15)，式(20)作为边界条件，通过求解以下最优化问题得到β₀₂.

(21)

算法描述见算法1.
算法1??迭代初值的选取
输入x，f_d1，k.
1) 数据x归一化；计算分类中心; X
2) 区间[-f_d1, f_d1]内随机生成β₀₁;
3) 以式(19)计算区间，在该区间内随机生成α₀;
4) 将α₀降序处理，判断式(16)是否成立.若是，转步骤5)；若否，转步骤2)；
5) 以式(15)，式(20)为约束，判断式(21)是否存在最优解.若是，得到β₀₂，结束；若否，转步骤2).
输出b₀=[α₀^T, β₀₁^T, β₀₂^T]^T.
通过此算法得到的b₀即可作为Newton法的迭代初值.
1.3 迭代过程的控制虽然得到了较好的迭代初值，但是在迭代过程中可能会出现如下情况：
1) α_k+1元素之间大小关系错乱；
2) β_k+1与β_k之间差距过大.
这些问题会导致Hessian矩阵奇异，迭代失败.通过分析，这些情况通常是由步长Δ_k= H^-1(b_k)·f(b_k) =[Δα_k; Δβ_k]过大造成.下面给出一种自适应控制方法.
当出现情况1)时，利用下式减小步长：

(22)

其中，w为调整系数，初值为1，一旦常系数在k+1次迭代不满足次序条件时便将w赋值为w/2，重新计算b_k+1.
为防止情况2)的出现，需重点控制Δβ_k+1元素之间的差距.设置阈值t_d，令u=max(|Δβ_k+1|).若u大于阈值t_d，则需要对其进行约束：

(23)

其中:放大系数f_d3推荐为1.2~1.5；阈值t_d为1~2.
算法描述见算法2.
算法2??迭代过程的控制
输入??b₀，f_d3，t_d，收敛阈值ε.
1) 初始化参数k=0，w=1;
2) 利用迭代格式(6)计算b_k+1;
3) 判断α_k+1元素间大小关系是否正确.若是，转步骤4)；若否，w= w/2，通过式(22)调整，转步骤3);
4) 判断u=max(|Δβ_k+1|)是否在阈值t_d之内.若是，转步骤5)；若否，通过式(23)调整，转步骤3);
5) ‖b_k+1 -b_k‖₂²是否小于ε.若是，结束；若否，k=k+1，w=1，转步骤2).
输出??b_k+1.
算法2中的“=”为赋值符号.由于算法1已保证了α₀的顺序，因此只要w调整适当，总会保持α_k的大小关系不变.
2 滚动轴承健康状态评估实例文献[9]通过非线性时间序列Logistic模型，验证了检测滚动轴承非线性突变信号的敏感性和有效性，因此本文也以轴承健康状态评估为例.选取美国凯斯西储大学轴承数据库^[5]中采样频率为12 kHz，负载为0，正常状态和滚动体故障(通过电火花加工设置故障，火花点直径为0.177 8，0.355 6，0.533 4 mm，分别模拟轻微磨损、一般磨损及严重磨损)的振动信号数据为原始数据，其中响应类别1，2，3，4分别表示正常状态、轻微磨损、一般磨损和严重磨损，模型建立过程如下.
首先，根据文献[10]，提取振动信号中的30个特征，得到相应特征值，再利用逐步模型选择法，以Wald统计量为指标从中逐步筛选出8个主要特征：F_C，F₂，F₅，F₆，F₇，F₈，F₉，F₁₂；然后，利用文献[11]所述剔除离群点方法，从所有606组数据中筛选出有效数据597组；最后，通过本文方法和SPSS软件分别建立累积Logistic模型，所得回归系数如表 1所示.
表 1(Table 1)

表 1 不同模型的回归系数Table 1 Regression coefficients of different models 模型 常系数 变系数 1 2 3 FC F2 F5 F6 F7 F8 F9 F12 本文 -3.915 -17.839 -31.791 -56.702 -3.054 84.526 8.477 -2 399.118 3 166.694 790.788 19.472 SPSS 0.781 -7.393 -14.495 -25.81 -2.565 39.929 6.188 -1 148.056 1 522.747 385.872 12.574

表 1 不同模型的回归系数 Table 1 Regression coefficients of different models

根据回归结果，两种方法所得各系数在组内所占权重大致相同，且各自变量系数的正负符号相同，说明两种方法对于当前数据具有相同的解释性.模型评价指标数值结果见表 2.
表 2(Table 2)

表 2 不同模型的评价指标Table 2 Evaluation indicators of different models 评价指标 本文模型 SPSS AIC 85.898 3 90.644 0 McFadden-R2 0.946 9 0.922 9 预测准确率 0.974 9 0.966 5 Gamma 0.999 6 0.999 2 Somers’D 0.986 4 0.982 4 Kappa 0.965 0 0.955 6

表 2 不同模型的评价指标 Table 2 Evaluation indicators of different models

表 2中除赤池信息量AIC外，其余指标的数值越大越能证明模型的优势.可以看出，本文方法在各项指标上都更加突出，证明了该方法的有效性.为进一步验证稳健性，分别利用两种方法进行200次Booststrap随机试验.每次试验在正常、轻微磨损、一般磨损及严重磨损这4种类别中，分别随机选取25组共100组数据作为验证样本，其余作为训练样本.将试验结果制成柱状分布图.两方法的训练准确率及验证准确率如图 1~图 4所示.
图 1(Fig. 1)

图 1 本文方法训练准确率分布图Fig.1 Distribution of training accuracy of the proposed method

图 2(Fig. 2)

图 2 SPSS训练准确率分布图Fig.2 Distribution of training accuracy of SPSS

图 3(Fig. 3)

图 3 本文方法验证准确率分布图Fig.3 Distribution of verification accuracy of the proposed method

图 4(Fig. 4)

图 4 SPSS验证准确率分布图Fig.4 Distribution of verification accuracy of SPSS

由图 1~图 4可知，两种方法训练准确率和验证准确率的分布都大致符合正态分布；通过对比可以直观看出，本文方法无论训练准确率还是验证准确率，都处于更高水平，表明该方法具有较强的稳健性.
3 结论1) 通过自适应地选取迭代初值和控制迭代过程，避免了Hessian矩阵奇异的情况，使Newton迭代法能够有效地进行累积Logistic模型的参数估计.
2) 与SPSS累积Logistic回归模型数值结果对比，本文方法在各项模型评价指标上具有更高的水平，验证了该方法的有效性.
3) 基于200次的Booststrap随机试验对比，表明本文方法对于数据的依赖性较小，具有较强的稳健性.
参考文献

[1]	Nizami N, Prasad N. Multinomial Logistic regression analysis[M]. Singapore: Palgrave Macmillan, 2017: 315-319.
[2]	Kanbayashi Y, Ishikawa T, Kanazawa M, et al. Predictive factors in patients eligible for pegfilgrastim prophylaxis focusing on RDI using ordered logistic regression analysis[J]. Medical Oncology, 2018, 33(5): 55-61.
[3]	Yu J B. Aircraft engine health prognostics based on Logistic regression with penalization regularization and state-space-based degradation framework[J]. Aerospace Science & Technology, 2017, 68: 345-361.
[4]	Boateng K A.Modeling external capital inflows and economic growth in Africa utilizing ordinal Logistic regression[D].Minneapolis: Walden University, 2013. http://search.proquest.com/docview/1468949031
[5]	Case Western Reserve University.Bearing data center[EB/OL].(2013-07-15)[2018-12-15].http://csegroups.case.edu/bearingdatacenter/home.
[6]	王济川, 郭志刚. Logistic回归模型——方法与应用[M]. 北京: 高等教育出版社, 2001: 237-242. (Wang Ji-chuan, Guo Zhi-gang. Logistic regression model—method and application[M]. Beijing: Higher Education Press, 2001: 237-242.)
[7]	Wang H Y.More efficient estimation for logistic regression with optimal subsample[EB/OL].(2018-03-31)[2018-12-15].https://arxiv.org/pdf/1802.02698.pdf.
[8]	Chiang L H, Kotanchek M E, Kordon A K. Fault diagnosis based on Fisher discriminant analysis and support vector machines[J]. Computers & Chemical Engineering, 2004, 28(8): 1389-1401.
[9]	刘永斌.基于非线性信号分析的滚动轴承状态监测诊断研究[D].合肥: 中国科学技术大学, 2011. (Liu Yong-bin.The diagnosis of rolling bearing condition monitoring based on nonlinear signal analysis[D].Hefei: University of Science and Technology of China, 2011. )
[10]	Lei Y, He Z, Zi Y, et al. Fault diagnosis of rotating machinery based on multiple ANFIS combination with GAS[J]. Mechanical Systems & Signal Processing, 2007, 21(5): 2280-2294.
[11]	Yamada M, Liu S, Kaski S.Interpreting outliers: localized Logistic regression for density ratio estimation[EB/OL].(2017-2-21)[2018-12-15].https://arxiv.org/pdf/1702.06354.pdf.