传统EM算法采用随机选取或经验选取的方法初始化待估计参数,容易造成迭代算法陷入局部最优。相关文献中采用绝对分类的初始化方法,如Kmeans分类、Binning分类。混合模型实际分布一般为多个独立分布的叠加,各独立分布之间存在交集,采用绝对分类初始化与理论情况不一致。本文提出一种贝叶斯随机分类的方法将不完全数据集随机分配到各独立分布,作为EM算法的初始化方法。贝叶斯随机分类方法步骤如下:

Step1 采用Kmeans聚类分析法对不完整数据集进行分类,对每类数据进行威布尔分布参数的极大似然估计,作为各独立分布的先验分布参数,第 j类样本数据对应的分布参数为( η_j, β_j)。每类样本数据容量 g_i占总容量 G的比值为各独立分布的权值 w_i,即 w_i=g_i/G。

Step2 计算第 i点在第 j类独立分布的密度函数值 f_ij( x_i),并根据式(6)归一化处理得到 p_ij, p_ij表示第 i个点隶属于第 j个独立分布的后验概率:

Step3 将第 i点分配到第 j类的后验概率区间为 ,其中p _i0=0。每个点按照均匀分布产生随机概率r _i,根据r _i是否属于第j类后验概率区间判别第i点所属的分类。并根据当前样本分类更新每个独立分布相应的参数。

Step4 根据式(7)计算最大似然函数值,选似然函数值较大者为当前最优分类:

Step5 重复Step2~Step4,直到最大似然函数值收敛或迭代次数结束。

EM算法的核心思想是:根据已有的数据,借助隐藏变量,通过期望值之间的迭代,估计似然函数。具体EM算法步骤如下:

Step1 初始化参数 , ,…, ; , ,…, ; , ,…, 。

Step2 期望步。根据式(8)迭代计算响应度 γ_i:

Step3 极大化步。根据方程式(9)迭代计算加权均值和方差:

Step4 重复步骤Step2、Step3直到极大似然函数值收敛。

在Step3中, 、 与 、 存在关系方程式(10),该方程显然为超越方程,难以得到解析解的形式。

在EM算法中,每一次迭代都要求解方程(10),令 λ=μ/σ,当 γ=0时,则有公式:

λ与 β显然存在单调关系,由于函数关系复杂,采用径向基函数插值的数值方法^{[ 13]}根据 λ求解 β。然后,根据式(12)求解 η:

为验证参数估计方法的准确性,采用逆变换分布函数的方法产生多重混合威布尔分布的样本数据。根据逆变换得 x=F^-1( y),其中 y服从(0,1)上的均匀分布。产生一组随机数 Y=[ y₁, y₂,…, y_n],根据逆变换表达式得到服从多重混合威布尔分布的样本数据 X=[ x₁, x₂,…, x_n]。给出不同重数的威布尔混合模型算例,分别产生容量大小为1000的随机样本,对比多重威布尔混合模型参数估计的常用方法,如传统EM算法、优化算法,给出实际算例的计算结果,如表1所示。

对比不同重数的混合威布尔分布算例,对比极大似然函数值 L_max,可以看出,在威布尔重数为2的情况下,各种算法与理论分布都比较接近,估计结果无明显差异。威布尔重数大于3的情况下,改进的EM算法的参数估计结果更接近理论分布。

表1

Table 1

表1(Table 1)

表1 不同算法对比参数估计值表 Table 1 Estimated parameters obtained by different methods 算例 2重威布尔分布算例 3重威布尔分布算例 4重威布尔分布算例 理论 分布 参数 改进 EM 算法 传统 EM 算法 模拟 退火 算法

表1 不同算法对比参数估计值表 Table 1 Estimated parameters obtained by different methods

威布尔混合模型分布密度函数形状多样,当形状参数为1时,代表指数分布;当形状参数为2时,代表瑞利分布;当形状参数为3时,与正态分布极为接近。可以考虑用多重威布尔混合分布逼近其他分布,如三参数威布尔分布、 Γ分布、对数正态分布等。建立多重威布尔分布逼近其他分布的最小二乘拟合模型:

式中: w_jf_j(x)表示 N重威布尔分布密度函数; g( x)表示待拟合其他分布密度函数。

采用4重威布尔分布分别逼近三参数威布尔分布、 Γ分布、对数正态分布、 χ²分布。逼近效果如图1所示,4重威布尔分布与所逼近的分布几乎完全重合,说明多重威布尔混合模型具有良好的逼近性能。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 多重威布尔混合模型逼近效果Fig.1 Approximation results of finite Weibullmixture distributions

采用传统的极大似然估计方法计算得到单重威布尔分布参数估计值,采用改进EM算法计算得到多重威布尔混合模型的参数估计值。参数取值如表3所示。

根据图2中不同重数的威布尔分布密度函数曲线,可以看出,随着威布尔混合模型重数的增加,数据分布规律趋于稳定。其中,4重与5重威布尔混合模型几乎重合,说明案例中的冲压机床可靠性分布规律在选用4重威布尔混合模型时达到稳定,更能反映该型号冲压机床早期可靠性的实际分布规律。

表3

Table 3

表3(Table 3)

表3 多重威布尔混合模型参数表 Table 3 Estimated parameters of finite Weibull mixture distributions 重数 N 权重 wi 尺度参数 ηi 形状参数 βi 1 w1 =1 η1 =363 .5 β1 =1 .64 2 w1 =0 .4835, w2 =0 .5165 η1 =522 .75, η2 =193 .905 β1 =4 .85, β2 =2 .48 3 4 5

表3 多重威布尔混合模型参数表 Table 3 Estimated parameters of finite Weibull mixture distributions

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 威布尔混合模型分布密度对比Fig.2 Density comparison of finite Weibullmixture distributions

拟合优度检验是检验来自总体中的一类数据的分布是否与某种理论分布相一致的统计方法。常用的检验法有 t检验、 χ²检验、KS检验。其中, t检验属于正态性检验方法,不适用于其他分布的检验。 χ²检验与KS检验相比,KS检验具有适用于小样本、检验结果的鲁棒性强等特点。因此,本文选用KS检验方法。

若 S_n( x)表示一个 n次观察的随机样本观察值的累积概率分布函数, S_n( x) =i/n, i是等于或小于 x的所有观察结果的数目, i=1,2,…, n。 F₀( x)表示一个特定的累积概率分布函数,也就是说,对于任一 x值, F₀( x)代表小于或等于 x值预期结果所占的比例。于是定义 S_n( x)与 F₀( x)之差的绝对值为 D,即:

式中: S_n( x)为经验分布函数; F₀( x)为理论分布函数,若对每一个 x值来说, S_n( x)与 F₀( x)差异很小,则表明经验分布函数与理论分布函数的拟合程度很高,认为样本数据来自具有该理论分布的总体。KS检验集中考察统计量 D_max,即:

KS检验步骤中,首先建立假设:

然后,计算统计量 D_max。根据给定的显著性水平 α,样本个数 n,查单样本KS检验统计量表得临界值 d_α。若 D_max< d_α,则在 α水平上,接受 H₀;否则,拒绝 H₀。根据该检验法得到混合威布尔分布模型重数 N对应的检验统计量 D_max为0 .0756( N=1)、0 .0457( N=2)、0 .0225( N=3)、0 .0278( N=4)、0 .0276( N=5)。

查单样本KS检验统计量表得显著性水平 α=0.05对应的临界值 d_0.05=0.2270。根据 D_max值,在显著性水平 α=0.05时,威布尔混合模型重数为1到5的情况下,原假设 H₀均被接受。并且,随着威布尔混合模型重数的递增,检验统计量 D_max呈递减趋势,说明威布尔混合模型重数越高,分布模型的拟合优度越好。