凸约束广义线性模型参数MLE估计的收敛性研究

钟绍军，童恒庆

（1.咸宁学院 数学与统计学院，湖北 咸宁 437100；2.武汉理工大学 数学与统计学院，武汉 430070）

1 凸约束广义线性回归模型及其参数估计

考虑以下的评估工作：假定有m个被评估对象，有p个评估指标，对每一个评估对象的每一个指标，都有n个专家进行评估打分。p个指标分别以x(1),…x(p)表示。第i(i=1,…,n)个专家对第k(k=1,…,m)个对象的第j(j=1,…p)个指标进行评估，得到的分数记为xijk，每个专家对每个对象的评估分数是一个p维向量数据，作为评估数据表的一行。显然，每个对象的评估数据有n行，构成一个数据块。m个数据块构成一个立体数据阵X={xijk}。给每个评估对象最终的评估分数只有一个，记为yk(k=1,…,m)，它也是待定的。评估模型可表述为

其中Dmn×n=Im⊗In，⊗是Kronecher积。这是童恒庆教授（1993）针对评估过程提出的一种特殊的凸约束广义线性回归模型[1]。童教授给出了模型参数的基于交互投影的最小二乘估计和基于EM算法的极大似然估计[2～4]。前者的收敛性得到了证明，但后者至今没有结果。本文主要解决模型（1）参数基于EM算法的极大似然估计的收敛性。

2 基于EM算法的参数极大似然估计[2][3]

假定模型（1）中关于随机误差ε～N(0,σ2I)，则得到对数似然函数为

一般假定σ2已知，而y是未知的。EM算法的E步即在（2）两边取关于y的条件数学期望，其密度函数f(y|β(i))（β(i)是β在第i次迭代的估计值）可以由随机误差的分布导出，于是得到

EM算法的M步即为求Q(β,β(i))的极大值。将条件密度函数代入化简，（3）的极大化相当于对下式的极小化：

记使q(β,β(i))在约束条件下取极小值的β的解为β(i+1)，获得β(i+1)就实现了EM算法的M步，将β(i+1)作为下一次迭代的初值代入（4）式，重复迭代即可得到序列，直到求解的误差小于预定精度＜ε)为止。

3 基于EM算法的极大似然估计的收敛性

若令A=(H'⋮X0')，b=(1,0,0,…,0)'，λ2=(k1,k2,…,kp)，λ=(λ1,λ2)'，λ1和ki(i=1,2,…,p)都是实数，则（6）式可转化为其对偶问题

其中，Aλ-g=θ，λ2≥0。若令A'θ-v=b，其中v为松弛变量，则

引理1对任意满足A'θ≥b的θ都有：

证明：由（7）式可得：A'λ(i)-g(i)=θ(i+1)，v(i)=A'Aλ(i)-A'g(i)-b，其中(v(i),λ(i))=0。于是有(v(i),λ(i))=0=(Aλ(i)Aλ(i))-(A'g(i),λ(i))-b'λ(i)。由A'λ(i)-g(i)=θ(i+1)，得到

(Aλ(i),Aλ(i))-(A'g(i),λ(i))=(A'θ(i+1),λ(i))

于是有 (A'θ(i+1),λ(i))=b'λ(i)。因为λ(i)≥0 ，故对任意满足A'θ≥b的θ都有：(λ(i),A'θ)=b'λ(i)。因而 (A'θ(i+1)-A'θ,λ(i))≤0 ，即 (θ(i+1)-θ,Aλ(i))≤0 ，由A'λ(i)-g(i)=θ(i+1)，可得(θ(i+1)+g(i),θ(i+1)-θ)≤0。引理得证。

定理1线性不等式约束下的EM算法序列{θ(i+1)-θ(i)}是收敛的并且收敛到零。

证明：先证明极大对数似然函数序列{l(θ(i)|y)}是有界的。由模型的参数的约束条件可得：β是满足配方条件1'pβ=1的非负的向量，从而有上界。因变量y是自变量X的配方汇总因而也是有上界的。那么观测对数似然函数{l(θ|y)}关于θ是有上界的。

由于g=-E(y|β(i))=-E(y|Λθ(i))是条件期望，它是关于θ(i)和X0的函数，故可将其分解为g=-Z-Cθ(i)，其中Z和C是关于X0的函数。那么EM算法的E步具有如下形式：

其中B(i)是关于θ(i)和X0的函数。那么由引理1可得

||θ(i+1)-θ(i)||2=||Z+Cθ(i)-θ(i)||2-||Z+Cθ(i)-θ(i+1)||2-2(Z+Cθ(i)-θ(i+1),θ(i+1)-θ(i))[≥||Z+Cθ(i)-θ(i)||2-||Z+C θ(i)-θ(i+1)||2]

又因：

由于序列{l(θ(i)|y)}是单调[5]的，且有上界，故必收敛。因此

于是可得线性不等式约束EM算法序列{θ(i+1)-θ(i)}是收敛的，并且

引理2 记M(⋅)为从θ(i)到θ(i+1)的映射,即θ(i+1)=M(θ(i)),那么必然存在一个向量C,对任意的θ1、θ2满足A'θ1≥b,A'θ2≥b, 都有

证明：显然，映射M(⋅)满足Lipschitz条件，且连续。

从映射M(⋅)的定义可知M(θ1)和M(θ2)满足线性不等式约束，即A'M(θ1)≥b，A'M(θ2)≥b。由引理1可得

(M(θ1)-M(θ2),Cθ1+Z-M(θ1)+M(θ2)-Cθ2-Z)≥ 0

于是有：||C||⋅||θ1-θ2||≥ ||C(θ1-θ2)||≥ ||M(θ1)-M(θ2)||

定理2线性不等式约束下的EM算法序列{}θ(i)存在收敛的子列，且{}θ(i)中任意收敛的子列都收敛到对数似然函数l(θ|y)在线性不等式约束A'θ≥b下的最优解。

证明：首先由于EM算法序列{}θ(i)满足Hθ(i)=1且θ(i)＞0，故{}θ(i)是有界的。那么必然存在收敛的子列，不妨设{}θ(ik)是{}θ(i)序列的任一收敛子列，并令其收敛点为θ*，即由引理2可得从而得到θ*也满足约束A'θ≥b.

由定理1知，映射M(⋅)是一个连续映射，故θ*=M(θ*)，从而可得

这说明θ*是Q(θ|θ*)在约束A'θ≥b下的一个K-T点，θ*为它的一个最优解。又由于Q(θ|θ*)是关于θ的严格凹函数，故由Kuhn-Tucker条件知Q(θ|θ*)在约束A'θ≥b下的最优解是唯一的。因此θ*是Q(θ|θ*)在约束A'θ≥b下的唯一最优解。而由θ*的唯一性，可知序列{}θ(i)是的任何收敛子列都收敛到Q(θ|θ*)在约束A'θ≥b下的最优解。

若θ1=θ2，则序列收敛。若θ1≠θ2，即θ1＞θ2，则序列必定存在两个收敛子列，满足

[1]Tong Hengqing.Evaluation Model and Its Iterative Algorithm by Alter⁃natingProjection[J].MathematicalandComputerModelling,1993,18(8).

[2]Hengqing Tong,Shaojun Zhong,Tianzheng Liu,Yanfang,Deng.Biosta⁃tistics Algorithm:Evaluation Model with Convex Constraint and Its Pa⁃rameters Estimates[C].The 1st International Conference on Bioinfor⁃matics and Biomedical Engineering,2007.

[3]童恒庆，余超，赵旭杰.凸约束广义线性回归模型的参数估计及算法[J].应用数学,2008,21(4).

[4]Dempster,A.P.,Laird,N.M.,Rubm,D.B.Maximum Likelihood from Incomplete Data Via the EM Algorithm(with Discussion)[J].Journal of the Royal Statistical Society,Series B，1977，（39）.