不同类型数据下混合模型参数估计效果的对比研究

Comparative Study on Effects of Parameter Estimation of Mixture Models under Different Types of Data

摘要: 混合高斯模型在描述数据方面应用较多，但它易受离群点的影响，其参数的极大似然估计不是稳健估计。混合t-分布模型由于其重尾分布的特性，相对于混合高斯分布，在分析重尾数据上更具稳健性。文章首先研究一元混合t-分布模型，利用标准EM算法给出了该模型参数极大似然估计的迭代步骤，并分别在三类模拟数据下与混合高斯模型进行了对比分析，验证了该模型的有效性以及在拟合重尾数据上的优势。算法初始化采用k-means方法。

Abstract: The normal mixture model has more applications in describing data. But it is easily influenced by the outlier, and the maximum likelihood estimation of parameters is not robust estimation. T-distribution mixture model has better robustness than Gauss mixture model to analyze data with longer time than normal tails because of its heavy-tails. In this paper, we studied a univariate t mixture model primarily. Based on EM algorithm, we derived the iteration steps of maximum li-kelihood estimation of the model’s unknown parameters. Furthermore, we did a comparative analysis by three types of simulated data. Simulation study shows that this model has an advantage in fitting data with longer time than normal tails. The initial value is given by k-means method.

1. 引言

混合分布模型是一种基于概率和统计的建模工具，它提供了一个用简单结构模拟复杂密度的灵活有效的方法，从而受到了统计学界、模式识别、图像处理等诸多领域的广泛关注。它的基本策略是，把研究数据看作是从很大的单个或多个总体上抽取出的一部分，通过潜在的分布或密度函数来描述。混合高斯模型应用的较多，但通常我们收集到的很多数据并不是严格的服从正态分布，而是较明显的服从重尾分布。混合t-分布模型由于其具有较长的尾巴，可对重尾点和异常点有效的降低权值，因此，相对于混合高斯分布模型，可以获得较强的精度和稳健性。

对于模型参数的求解，1977年Dempster等人在文献 [1] 中提出的EM算法成为了混合模型参数极大似然估计的极有效的工具。Peel和McLachlan在文献 [2] 中指出EM算法可以获得有限混合任意分布模型的极大似然估计。对于单一的t-分布，为了使M步更好求解，Meng和Rubin在文献 [3] 中用受限制的最大化CM步来代替M步，得到ECM算法；Peel和McLachlan在文献 [2] [4] 中提出多元混合t-分布模型，用标准EM算法求解混合t-分布模型参数的极大似然估计，并给出了ECM算法的一个应用；而后，Liu和Rubin在文献 [5] 中对ECM算法进行两处修改，得到收敛速度更快的ECME算法。在算法初始化方面，冉延平在文献 [6] 中用k-means方法确定混合高斯分布的最大混合子分布数目以及混合比例；史鹏飞在文献 [7] 中通过k-means方法先给出混合数据的一个粗糙分组，然后根据分组数据给出参数的一个粗略估计值，作为混合高斯分布EM算法的迭代初始值。随着计算机性能的快速发展，混合t-分布模型已应用于诸多领域。如杨云飞在文献 [8] 中提出了自适应均值滤波的多元t-分布混合模型，对医学脑图像分割进行了研究；熊太松在文献 [9] 中对伯克利图像数据和微软剑桥研究院提供的图像集用视觉和量化对比的评估方式证明了基于空间平滑的t-分布混合模型在灰度图像分割中的有效性；朱志娥在文献 [10] 中针对偏t正态数据、异方差和线性回归提出了偏t正态数据下混合线性联合位置与尺度模型，详细介绍了该模型下的EM算法并进行了有效的模拟验证。

在前人研究的基础上，本文研究了基于EM算法的三总体一元混合t-分布模型参数的极大似然估计，克服了多元混合t-分布模型中协方差矩阵向一元混合t-分布模型中尺度参数转变的困难，并首次将k-means方法用于该模型下算法初值的选取。此外，我们引进了混合高斯模型，并分别在三种不同类型数据下进行了对比仿真实验，验证了本文研究的模型和方法的有效性及其在处理重尾数据上的优势。

2. 有限t-分布混合模型

2.1. 一元学生t-分布

设随机变量y服从一元学生t-分布，记做 $y\sim t\left(y|\mu ,\sigma ,\nu \right)$ ，概率密度函数定义为 [11] ：

$t\left(y|\mu ,\sigma ,\nu \right)=\frac{\Gamma \left(\frac{1+\nu }{2}\right)}{\Gamma \left(\frac{\nu }{2}\right)}{\left(\frac{1}{\text{π}{\sigma }^2\nu }\right)}^{\frac{1}{2}}{\left[1+\frac{{\left(y-\mu \right)}^2}{{\sigma }^2\nu }\right]}^{-\left(\frac{1+\nu }{2}\right)}$ (1)

其中参数 $\mu$ 和 $\sigma$ 分别表示t-分布的位置参数和尺度参数， $\Gamma (·)$ 表示伽马函数，参数 $\nu$ 称为t-分布的自由度。

2.2. t-分布有限混合模型

设随机向量服从一元t-分布且由m个子分布混合而成，混合比例为 ${\pi }_k$ 且满足 ${\displaystyle \underset{k=1}{\overset{m}{\sum }}{\pi }_k}=1$ ， $t\left(y|{\theta }_k\right)$ 为第

k个子分布的概率密度函数， ${\theta }_k=\left\{{\mu }_k,{\sigma }_k,{\nu }_k\right\}$ 为对应的子分布的概率密度函数参数， $\Psi =\left\{{\pi }_1,\cdots ,{\pi }_m,{\theta }_1,\cdots {\theta }_m\right\}$ 为参数空间。因此，有限混合t-分布模型可以表示为

$t\left(y|\theta \right)={\displaystyle \underset{k=1}{\overset{m}{\sum }}{\pi }_{k}t\left(y|{\theta }_k\right)}$ (2)

本文研究三个子分布的情况，即取 $m=3$ ，则式(2)化为

$t\left(y|\theta \right)={\pi }_{1}t\left(y|{\theta }_1\right)+{\pi }_{2}t\left(y|{\theta }_2\right)+\left(1-{\pi }_1-{\pi }_2\right)t\left(y|{\theta }_3\right)$ (3)

其中， $t\left(y|{\theta }_k\right)$ 为第k个子分布的概率密度函数，具体形式见式(1)。

3. 模型参数极大似然估计的EM算法

本文要研究的模型为上文所提到的式(3)，并假设三个子分布的自由度相同，即 ${\nu }_1={\nu }_2={\nu }_3=\nu$ 。我们采用标准EM算法来求解模型参数，它提供了一种近似计算含有隐变量概率模型的极大似然估计的方法。在EM算法的基本框架下，我们引入隐变量以得到完整数据集。完整数据集定义为 $Y_c=\left\{Y,Z,U\right\}$ ，

其中， $Z$ 为标签变量 $Z=\{z_1,z_2,\cdots ,z_N\}$ ，且 $z_{ik}=\{\begin{array}{l}1,第i个样本来自第k个子分布\\ 0,其他\end{array}$ ( $i=1,2,\cdots ,N$ , $k=1,2,3$ )；

$U$ 为引进的另一个隐变量 $U=\left\{u_1,u_2,\cdots ,u_N\right\}$ ，给定 $z_{ik}=1$ 时， $u_i$ 是独立同分布的，且满足

$u_i|z_{ik}=1\sim G\left(\frac{\nu }{2},\frac{\nu }{2}\right)$ ； $Y$ 为可观测数据集 $Y=\left\{y_1,\cdots ,y_N\right\}$ ，且有 $y_i|u_i,z_{ik}=1\sim N\left({\mu }_k,{\sigma }_k^2/u_i\right)$ 。由(1) (3)两式

建立完整数据的似然函数为：

$\begin{array}{c}p\left(Y,Z,U|\Psi \right)={\displaystyle \underset{i=1}{\overset{N}{\prod }}p\left(y_i,z_{i1},z_{i2},u_i|\Psi \right)}\\ ={\displaystyle \underset{i=1}{\overset{N}{\prod }}{\displaystyle \underset{k=1}{\overset{3}{\prod }}{\left\{\frac{\sqrt{u_i}}{\sqrt{2\text{π}}{\sigma }_k}{\text{e}}^{-\frac{u_i{\left(y_i-{\mu }_k\right)}^2}{2{\sigma }_k^2}}\left[\frac{{\left(\frac{\nu }{2}\right)}^{\frac{\nu }{2}}}{\Gamma \left(\frac{\nu }{2}\right)}{\text{e}}^{-\frac{\nu }{2}{u}_i}{u}_i{}^{\frac{\nu }{2}-1}\right]{\text{π}}_k\right\}}^{z_{ik}}}}\end{array}$ (4)

则完整数据的对数似然函数为：

$\begin{array}{l}\log p\left(Y,Z,U|\psi \right)={\displaystyle \underset{k=1}{\overset{3}{\sum }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{z}_{ik}\{\frac{1}{2}\log u_i-\frac{1}{2}\log 2\text{π}-\log {\sigma }_k-\frac{u_i}{2{\sigma }_k^2}{\left(y_i-{\mu }_k\right)}^2}}\\ +\frac{\nu }{2}\log \frac{\nu }{2}-\log \Gamma \left(\frac{\nu }{2}\right)-\frac{\nu }{2}{u}_i+\log {\pi }_k\}\\ ={\displaystyle \underset{k=1}{\overset{3}{\sum }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{z}_{ik}\log {\pi }_k}}\\ +{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{z}_{ik}\left\{-\log \Gamma \left(\frac{\nu }{2}\right)+\frac{\nu }{2}\log \left(\frac{\nu }{2}\right)+\frac{\nu }{2}\left(\log u_i-u_i\right)\right\}}}\\ +{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{z}_{ik}\left\{-\frac{1}{2}\log 2\text{π}-\frac{1}{2}\log u_i-\log {\sigma }_k-\frac{u_i{\left(y_i-{\mu }_k\right)}^2}{2{\sigma }_k^2}\right\}}}\end{array}$ (5)

EM算法是一种迭代近似求解算法，它主要分两步进行：E步是对对数似然函数求期望，M步是最大化对数似然函数以获得新的参数值。

应用EM算法于上式，求解第 $j$ 次各参数的极大似然更新表达式。

E步： $Q\left(\Psi |{\Psi }^{\left(j\right)}\right)=E_{Z,U}\left[\log p\left(Y,Z,U|\Psi \right)|Y,{\Psi }^{\left(j\right)}\right]$

首先计算关于隐变量 $Z$ ， $U$ 的条件分布的期望：

$E_{Z,U}\left[z_{ik}|Y,{\Psi }^{\left(j\right)}\right]={\tau }_{ik}^{\left(j\right)}=\frac{{\pi }_k^{\left(j\right)}{f}_k\left(y_i|{\theta }_k^{\left(j\right)}\right)}{{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\pi }_k^{\left(j\right)}{f}_k\left(y_i|{\theta }_k^{\left(j\right)}\right)}}$ (6)

$E_{Z,U}\left[u_{ik}|y,{\Psi }^{\left(j\right)}\right]=u_{ik}^{\left(j\right)}=\frac{{\nu }^{\left(j\right)}+1}{{\nu }^{\left(j\right)}+{\left(y_i-{\mu }_k^{\left(j\right)}\right)}^2/{\sigma }_k^{2\left(j\right)}}$ (7)

$E_{Z,U}\left[\log u_{ik}|y,{\Psi }^{\left(j\right)}\right]=l_{ik}^{\left(j\right)}=\log u_{ik}^{\left(j\right)}+\left\{\psi \left(\frac{{\nu }^{\left(j\right)}+1}{2}\right)-\log \left(\frac{{\nu }^{\left(j\right)}+1}{2}\right)\right\}$ (8)

其中， $\psi \left(x\right)=\frac{\text{d}\left(\log \Gamma \left(x\right)\right)}{\text{d}x}=\frac{{\Gamma }^{\prime }\left(x\right)}{\Gamma \left(x\right)}$ ，于是

$\begin{array}{l}Q\left(\Psi |{\Psi }^{\left(j\right)}\right)={\displaystyle \underset{k=1}{\overset{3}{\sum }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}\log {\pi }_k}}\\ +{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}\left\{-\log \Gamma \left(\frac{\nu }{2}\right)+\frac{\nu }{2}\log \left(\frac{\nu }{2}\right)+\frac{\nu }{2}\left(\log u_{ik}^{\left(j\right)}+\Psi \left(\frac{{\nu }^{\left(j\right)}+1}{2}\right)-\log \left(\frac{{\nu }^{\left(j\right)}+1}{2}\right)-u_{ik}^{\left(j\right)}\right)\right\}}}\\ +{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}\left\{-\frac{1}{2}\log 2\text{π}-\log {\sigma }_k-\frac{u_{ik}^{\left(j\right)}{\left(y_i-{\mu }_k\right)}^2}{2{\sigma }_k^2}-\frac{1}{2}\left[\log u_{ik}^{\left(j\right)}+\Psi \left(\frac{{\nu }^{\left(j\right)}+1}{2}\right)-\log \left(\frac{{\nu }^{\left(j\right)}+1}{2}\right)\right]\right\}}}\end{array}$ (9)

M步： ${\theta }^{\left(j+1\right)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(j))}$

利用Q函数对各参数求偏导数并令其等于零，求解得到各参数的第 $j+1$ 次迭代更新表达式：

${\pi }_k^{\left(j+1\right)}={\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}/N}$ (10)

${\mu }_k^{\left(j+1\right)}=\frac{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}{u}_{ik}^{\left(j\right)}{y}_i}}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}{u}_{ik}^{\left(j\right)}}}$ (11)

${\sigma }_k^{\left(j+1\right)}=\text{sqrt}\left(\frac{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}{u}_{ik}^{\left(j\right)}{\left(y_i-{\mu }_k^{\left(j+1\right)}\right)}^2}}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}{u}_{ik}^{\left(j\right)}}}\right)$ (12)

自由度 ${\nu }^{\left(j+1\right)}$ 是非线性方程

$\log \frac{\nu }{2}-\psi \left(\frac{\nu }{2}\right)+1+\frac{1}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}}}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\tau }_{ik}^{\left(j\right)}\left(l_{ik}^{\left(j\right)}-u_{ik}^{\left(j\right)}\right)}=0$ (13)

的解。上式是关于 $\nu$ 的非线性方程，文献 [5] 中采用搜索 $\nu$ 的空间求出 $\nu$ 的估计值，但计算量大。文献 [12] 中作者给出了一个计算量相对较小的可直接计算多元混合t-分布模型中参数 $\nu$ 近似解的方法，但其不稳定。因此，方便起见，在接下来的数值模拟中我们只考虑自由度是已知的情形，即每一子分布的自由度都相同且固定为 $\nu$ 。

4. 数值模拟

为了验证上述参数估计方法的有效性，我们共采用三种不同类型数据进行模拟研究，并引进混合高

斯分布模型 [13] 与之作对比。由t-分布的方差与尺度参数的关系 ${\upsilon }^2=\frac{\nu }{\nu -2}{\sigma }^2$ ，将混合t-分布EM算法参

数估计结果中的尺度参数 $\sigma$ 转化为标准差 $\upsilon$ ，再与混合高斯分布EM算法估计的参数 $\upsilon$ 作比较。算法的初始化均采用k-means方法。参数估计的精确度采用均方误差来衡量，如混合比例 ${\pi }_1$ 的均方误差定义为：

$MSE\left({\pi }_1\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left({\pi }_{1i}-{\pi }_{1\left(0\right)}\right)}^2}$

其中， ${\pi }_{1\left(0\right)}$ 是 ${\pi }_1$ 的真值，n为模拟次数。

4.1. 混合高斯分布数据

给定真值 ${\pi }_{1\left(0\right)}=0.5$ 和 ${\pi }_{2\left(0\right)}=0.3$ ， ${\mu }_{1\left(0\right)}=2$ 、 ${\mu }_{2\left(0\right)}=7$ 、 $$ 、 ${\upsilon }_{1\left(0\right)}={\upsilon }_{2\left(0\right)}={\upsilon }_{3\left(0\right)}=1$ ，分别取样本量 $N=500$ ，1000，共产生2组混合高斯分布数据。对混合t-分布模型，分别取自由度 $\nu =3$ [14] ，15，30。重复模拟1000次，模拟结果如下：

由表1、表2、表3知：

$\nu =3$ 时，混合高斯模型参数估计的均方误差均比混合t-分布模型参数估计的均方误差小，这一点在 ${\upsilon }_1$ 、 ${\upsilon }_2$ 、 ${\upsilon }_3$ 上更为明显；在 $\nu =15,30$ 时，除了对 ${\upsilon }_1$ 、 ${\upsilon }_2$ 、 ${\upsilon }_3$ 的估计结果混合高斯模型略好于混合t-分布模型外，两种方法对其他参数的估计的均方误差，几乎无差。此外还有，随着自由度的增大，混合t-分布模型参数估计的均方误差变小；整体来看，样本量越大，MSE越小。

4.2. 混合t-分布数据

给定真值 ${\pi }_{1\left(0\right)}=0.5$ 和 ${\pi }_{2\left(0\right)}=0.3$ ， ${\mu }_{1\left(0\right)}=2$ 、 ${\mu }_{2\left(0\right)}=7$ 、 ${\mu }_{3\left(0\right)}=11$ 、 ${\upsilon }_{1\left(0\right)}={\upsilon }_{2\left(0\right)}={\upsilon }_{3\left(0\right)}=1$ ，取 $\nu =3$ [14] ，15，30。分别取样本量 $N=500$ ，1000，共产生6组混合t-分布数据。重复模拟1000次，模拟结果见下表：

由表4、表5、表6知：

混合t-分布模型可以较好地拟合该数据，参数估计值与真值十分接近。当 $\nu =3$ 时，对所有参数的估计，混合t-分布模型参数估计的均方误差均比混合高斯分布模型参数估计的均方误差小； $\nu =15$ 时，除 ${\mu }_2$ 、

Table 1. The simulation results of ν = 3

表1. $\nu =3$ 的模拟结果

Table 2. The simulation results of ν = 15

表2. $\nu =15$ 的模拟结果

Table 3. The simulation results of ν = 30

表3. $\nu =30$ 的模拟结果

注：EST1：混合高斯分布模型下的平均估计值。EST2：混合t-分布模型下的平均估计值。MSE1：混合高斯分布模型参数估计的均方误差。MSE2：混合t-分布模型参数估计的均方误差。

Table 4. The simulation results of ν = 3

表4. $\nu =3$ 的模拟结果

Table 5. The simulation results of ν = 15

表5. $\nu =15$ 的模拟结果

Table 6. The simulation results of ν = 30

表6. $\nu =30$ 的模拟结果

${\upsilon }_2$ 、 ${\upsilon }_3$ 外，混合t-分布模型参数估计的均方误差比混合高斯分布模型参数估计的均方误差小；在 $\nu =30$ 时，除 ${\upsilon }_3$ 外，混合t-分布模型参数估计的均方误差均比混合高斯分布模型参数估计的均方误差小。此外，随着自由度的增大，混合t-分布模型参数估计的均方误差变小；整体来看，样本量越大，MSE越小，估计结果越好。

4.3. 含噪声的混合高斯数据

混合t-分布模型相对于混合高斯模型有着较好的稳健性，这种稳健性尤其体现在对重尾数据(含噪声点、异常点数据)的处理。而处理重尾数据的另一种方法是在高斯分布的基础上添加一个均匀分布的成分 [6] 。因此，我们在双高斯数据的基础上添加一个均匀分布的部分作为重尾数据，再分别用混合t-分布模型和混合高斯模型进行拟合并作比较。因为前两小节已经对自由度、样本量进行了研究比较，在这一小节我们不再考虑此二者的影响。取噪声所占比例分别为5%和10%，混合比例 ${\pi }_1=0.5$ 、 ${\pi }_2=0.3$ ，自由度 $\nu =15$ ，样本量 $N=1000$ 。重复模拟1000次，模拟结果如表7和表8。

由表7、表8知：

通过比较两种模型下参数的估计结果和均方误差我们可以得到，混合t-分布模型对该类型数据拟合的较好，尤其对混合比例、位置参数的估计都较混合高斯分布模型估计的效果好。而对于尺度参数的估计，混合高斯模型拟合下得到的参数的均方误差略小，但相差不大。因此相对于混合高斯分布，混合t-分布模型可以更好的拟合含噪声的混合高斯数据，这也正说明了混合t-分布模型较于混合高斯模型能够更好地处理重尾数据。

Table 7. Mixed Gaussian data with 5% noises

表7. 含噪声5%的混合高斯数据

Table 8. Mixed Gaussian data with 10% noises

表8. 含噪声10%的混合高斯数据

5. 结论

本文主要研究了一元混合t-分布模型，给出了EM算法下该模型参数的极大似然估计，并在模拟的三种类型的数据下与混合高斯模型进行了对比分析。从前两类数据的模型参数估计结果中可以看出，每个子分布的自由度固定且取相同的值的情况下，对于混合高斯数据，当自由度的取值足够大时，基于混合t-分布模型的EM算法的参数估计结果并不比基于混合高斯模型的EM算法差；对于混合t-分布数据，基于混合t-分布模型的EM算法能够得到较好的估计结果并优于基于混合高斯模型的EM算法的估计结果，且随着自由度的增大，效果会更好；而在第三类含噪声的混合高斯分布数据下，整体而言，混合t-分布模型比混合高斯分布模型拟合效果更好，说明了混合t-分布模型在处理重尾数据上更具优势。以上结果验证了本文研究的模型和方法的有效性。

资助信息

国家自然科学基金(11601150)；国家自然科学基金(U1430103)；中央高校基本科研业务费专项资金资助(2016MS62)。