Application of EM Algorithm in Parameter Estimation of p-Norm Mixture Model

设$l_i(i=\mathrm{1,\; 2}\cdots n)$是$\mu$的一组独立观测值，观测值$L=(l_{\mathrm{1}},l_{\mathrm{2}}\cdots l_n{)}^{\mathrm{T}}$服从一元p范分布，概率密度函数为^[1，3]：

式中，$\lambda =\sqrt[]{\Gamma (\mathrm{3}/p)/\Gamma (\mathrm{1}/p)}$；参数分量$\theta =(p,\mu ,\sigma )$，其中，p为尺度参数（当p=1时，误差的概率密度函数为拉普拉斯分布；当p=2时，误差的概率密度函数为正态分布；当p值无限趋近于0时，误差服从极限分布；当p值无限趋近于正无穷大时，误差服从均匀分布）。

若各观测值服从一元p范分布，假定模型的混合数为$m$，则构成了一元p范混合模型，该模型的概率密度函数表示如下：

式中，参数分量$\theta =({\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}}\cdots {\theta }_i\cdots {\theta }_m)$；${\theta }_i=(p_i,{\mu }_i,{\sigma }_i)$；${\alpha }_j$为混合的权重，表示满足第$j$种分布的数据所占的比例，为了满足密度函数的性质，必须满足$\mathrm{0}\le {\alpha }_j\le \mathrm{1}$且${\displaystyle \sum _{j=\mathrm{1}}^m}{\alpha }_j=\mathrm{1}$，第$j$个总体的密度函数为$f_j(L)$。

该混合模型具有较强的灵活性，其难点在于如何求解参数的估计量。p值可以根据一定的迭代方法进行确定^[1，3]，也可以采用直接计算公式快速估计出p值^[2，15]，从而待估计的参数只有方差${\sigma }_i$、均值${\mu }_i$以及混合数${\alpha }_i$。

设$L$为混合分布的观测数据，由于无法分辨出哪个样本来自哪个分布，因此，观测数据中没有包含数据的全部信息，是不完全数据。引入分量$z_{ij}$，当$z_{ij}=\mathrm{1}$时，表示第$i$个观测数据来自第$j$个分布，即$P(z_{ij}=\mathrm{1})={\alpha }_j$；当$z_{ij}=\mathrm{0}$时，表示第$i$个观测数据不是来自第$j$个分布，则观测值$L$的条件分布密度函数为：

由于混合数${\alpha }_j$是无法观测的，因此称为不完全数据或缺失数据，设缺失数据向量为$\mathit{Z}=(z_{\mathrm{1}},z_{\mathrm{2}}\cdots z_n{)}^{\mathrm{T}}$，$z_i=(z_{i\mathrm{1}},z_{i\mathrm{2}}\cdots z_{im})$，则$(L,\mathit{Z})$称为完整数据。则有：

完全数据$(L,\mathit{Z})$的似然函数为：

缺失数据$Z$的条件分布为：

设$l_{\mathrm{1}},l_{\mathrm{2}}\cdots l_n$为取自上述p范混合模型的一组独立观测数据，对应的缺失数据为$z_{\mathrm{1}},z_{\mathrm{2}}\cdots z_m$，记$\theta =(\alpha ,\mu ,\sigma )$，则完全数据的对数似然函数为：

对含有$m$个子体的p范混合模型来说，EM算法是迭代算法。先给定参数的初始值${\alpha }_{\mathrm{1}}^{(\mathrm{0})}$，${\alpha }_{\mathrm{2}}^{(\mathrm{0})}$$\cdots$${\alpha }_m^{(\mathrm{0})}$；$({\mu }_{\mathrm{1}}^{(\mathrm{0})},{\sigma }_{\mathrm{1}}^{(\mathrm{0})})$，$({\mu }_{\mathrm{2}}^{(\mathrm{0})},{\sigma }_{\mathrm{2}}^{(\mathrm{0})})$$\cdots$$({\mu }_m^{(\mathrm{0})},{\sigma }_m^{(\mathrm{0})})$，由它求出缺失数据的值；再根据此数据估计出新的参数估计值，根据这一估计值对缺失数据的值进行更新；如此反复迭代，直到收敛为止。

应用EM算法求解式（7）。求解第k次各参数表达式的步骤如下：

1）E步：构造$Q(\theta ,{\theta }^{(k)})$。

缺失数据$\mathit{Z}$的条件分布期望为：

由式（7）、（8）得：

式中，$W_{ij}^{(k)}=({\alpha }_j^{(k)}{f}_j^{(k)}(l_i))/(f^{(k)}(l_i))$。

2）M步：将对数似然函数$Q(\theta ,{\theta }^i)$极大化，求取相应参数的参数估计值。在${\displaystyle \sum _{j=\mathrm{1}}^m}{\alpha }_j=\mathrm{1}$的限制条件下，由式（7）求$Q(\theta ,{\theta }^{(k)})$关于${\alpha }_j$、${\mu }_j$和${\sigma }_j$的最大值，对似然函数各参数求导，令其等于0，有：

由于${\displaystyle \sum _{j=\mathrm{1}}^m}{\alpha }_j=\mathrm{1}$，可以得到$\lambda ={\displaystyle \sum _{j=\mathrm{1}}^m}{\displaystyle \sum _{i=\mathrm{1}}^n}{W}_{ij}=n$，故得到${\alpha }_j$的估计值，从而得到该参数的第$k+\mathrm{1}$次迭代表达式为：

分别对${\mu }_j$、${\sigma }_j$求偏导，化简，同理可得到参数的第$k+\mathrm{1}$次迭代更新式如下：

式（12）~（14）是参数的非线性方程，可以采用迭代的方法计算，计算步骤可总结为：（1）‍选择合适的初始值，令${\theta }_j^{(k)}=({\alpha }_j^{(k)},{\mu }_j^{(k)},{\sigma }_j^{(k)})$；（2）进行第k+1次迭代，求得新的混合系数${\alpha }_j^{(k+\mathrm{1})}$；（3）计算均值的估计值，通过迭代解方程‍（13），求得第k+1次迭代的均值${\mu }_j^{(k+\mathrm{1})}$；（4）计算方差的估计值，将第‍（2）步、第‍（3）步中得到的混合系数${\alpha }_j^{(k+\mathrm{1})}$和${\mu }_j^{(k+\mathrm{1})}$代入式（14），求得第k+1次迭代的方差值${\sigma }_j^{(k+\mathrm{1})}$；（5）比较迭代后得到的各参数估计值与迭代前相应参数估计值的差值是否充分小，若不满足，则将此次迭代值作为下一次迭代的初始值进行迭代运算，直到差值充分小停止循环。

假定模型的混合系数为2，采用拉普拉斯分布子样、高斯分布子样以及实测GPS观测值残差数据作为实验数据。首先，利用矩估计求解出混合数据的p值^[16]；然后，利用本文的EM算法解算混合p范分布的参数值；最后，通过分析实验的结果验证本文算法的可行性。

假定p范混合模型为高斯混合模型，即混合数据由两组高斯分布数据组成。通过Matlab软件随机生成服从$L_{\mathrm{1}}~N(\mathrm{3,\; 1})$，$L_{\mathrm{2}}~N(\mathrm{4},{\mathrm{3}}^{\mathrm{2}})$的数据，L₁中取600个数据作为样本1，L₂中取400个数据作为样本2，即样本总数$n=\mathrm{1}\mathrm{000}$进行实验。样本1和样本2的数据分布直方图如图 1所示（横坐标和纵坐标均无单位）。

Figure 1.  Sample Histogram and Distribution Curves of Example 1

利用矩估计迭代法求解出混合数据的模型分布参数p值，并以直方图的形式给出模型参数p值求解的可靠性，如图 2所示。

Figure 2.  Histogram and True and Estimated Distribution Curves (Example 1)

从图 2可以看出，两种高斯数据混合后，计算出的p值为1.363 9，可见混合高斯分布的数据不再服从高斯分布。与高斯分布相比，采用p范分布（p=1.363 9）估计出的概率密度与真实密度更为一致，也更加符合实际分布的情况。

为了获得更高的模型参数精度，利用上述样本数据进行12次实验（p选取1.363 9），将本文算法EM_p估计的结果列入表 1，取12次结果的平均值作为模型参数的最终值。为了说明混合模型参数估计的准确性，将计算结果与EM解算高斯混合模型参数（EM_G）结果进行对比，并采用均方根误差（root mean square error，RMSE）来表征参数估计的效果。

通过表 1可以发现，不论是EM_p还是EM_G算法，均能较好地估计出混合模型的6个参数。但从最终估计的结果来看，EM_p估计出的6个模型参数十分接近真值，估计精度远远高于EM_G算法。同时，每次估计出的模型参数变化均较小，且与真值符合度较高，从而验证了本文算法估计混合多峰数据的有效性和稳定性。结合图 3可以看出，EM_p算法估计的模型参数的RMSE均小于0.05，远远优于EM_G算法，进一步说明利用EM_p算法估计的混合高斯分布模型参数具有较好的精确度和稳定性。

Figure 3.  Root Mean Square Error of Model Parameter Estimates(Example 1)

假定混合模型中p值分别取1和2，混合数据由一组拉普拉斯分布数据$(L_{\mathrm{1}}~L(u=\mathrm{0},\sigma =\mathrm{1}))$和一组高斯分布数据$L_{\mathrm{2}}~N(\mathrm{1},{\mathrm{4}}^{\mathrm{2}})$组成。分别从L₁和L₂中各取1 000个数据作为样本数据，即样本总数$n=\mathrm{2}\mathrm{000}$进行实验，样本直方图见图 4。

Figure 4.  Sample Histogram and Distribution Curves of Example 2

利用矩估计法求解出混合数据的模型参数p值，由图 5可以看出，当p取0.952 6时，估计出的概率密度与真实密度十分接近，说明此时的模型更加符合该实验样本数据分布的真实情况。

Figure 5.  Histogram and True and Estimated Distribution Curves (Example 2)

表 2统计了混合数据的模型参数，取12次计算结果的平均值作为模型参数的最终值。

通过表 2可以发现，当样本数据是由拉普拉斯分布与高斯分布混合组成时，EM_G算法估计出的模型参数精度较差，不能将混合数据分类出来，其原因是样本数据中存在不同于高斯分布的数据，致使其算法失效。同时可以看出EM_p算法虽然模型参数估计的精度不如混合同分布数据时那么精确，但从最终估计的结果来看，EM_p算法仍能够较好地估计出混合模型的6个参数，估计出的模型参数十分接近真值，估计精度也远远高于EM_G算法。结合图 6可以看出，EM_G算法估计的模型参数的RMSE在4左右，估计的精度达不到模型参数估计所需的精度，而EM_p算法估计的模型参数的均方根误差在0.06附近，精度较高，充分说明EM_p算法能够有效地估计出混合异分布模型的参数。

Figure 6.  Root Mean Square Error of Model Parameter Estimates (Example 2)

数据来自加拿大Algonquin Park的ALGO测站点，利用TPS NET-G3A接收机采集获得2013-04-28的观测数据。在获得的32 ‍颗卫星对地观测数据中，选取某颗卫星伪距的精密单点定位双频无电离层组合观测值残差进行分析。取其中200 ‍个误差值作为样本数据，利用矩估计求出样本数据的p值为1.398。假设样本数据由两种分布数据组成，利用EM_p算法进行参数解算，结果如表 3所示。

由表 3可以看出，EM_p算法计算出的两类分布数据十分相似，几乎为同一种分布。所以将其模型分量数设为1，重新利用EM_p算法进行解算，取12个计算结果的平均值作为真值。最终求得GPS观测值残差服从p为1.398、μ为-‍0.052、σ为0.446的p范分布。

通过GPS伪距单点定位的精度来验证计算结果的正确性。利用武汉大学精密单点定位软件对观测数据进行解算，得到观测站点的精密坐标值，并以此作为伪距单点定位中测站点的真实坐标。现以p_0371180点的观测数据为例，对其进行解算，通过精密单点定位软件解算得到的测站点高精度坐标为：[X，Y，Z]=[-‍1 ‍304 ‍152.045 9，-4 831 831.378 5，3 ‍943 ‍232.966 1]。分别假设误差服从高斯分布和p范分布（p=1.398），利用最小二乘（least square，LS）和p范平差（least p-norm adjustment，Lp）求解伪距单点定位误差方程，得到待求参数的估计值。计算出定位点的三维坐标估计值，采用各历元解得的坐标估计值与真实坐标值进行对比，以每个方向上的坐标中误差作为GPS伪距单点定位的精度，见表 4。

从表 4可以看出，假定GPS观测值的误差服从高斯分布，利用传统的LS求解伪距单点定位的精度较低，其定位误差（1.138 m）达到米级。采用本文EM算法求解出的GPS观测值的误差分布模型进行p范平差所得到的坐标在3 ‍个方向上的精度均优于LS，定位效果达到分米级，从而进一步验证了本文EM_p算法估计模型参数具有较高的准确度与可靠性。

p范混合模型作为一种新的分布模型，考虑了测量误差的不确定性和误差分布的多样性。本文探讨了该模型参数的EM算法的迭代公式，利用模拟数据验证了EM算法结合p范分布可以有效解决误差的参数估计问题，并将其应用到GPS观测值误差分析中。从仿真和实测数据的算例计算结果可以看出，相比高斯混合模型，p范混合模型能够更好地反映出混合数据的实际分布情况，同时利用EM算法求解出的混合p范模型的参数值也更加准确。

本文扩充了p范分布理论，对进一步提高测量数据处理的精度具有一定的实用价值。同时还发现混合数的确定、尺度参数的确定和参数的初始值对计算结果影响较大，如何减小其影响是下一步需要考虑的问题。