考虑具有D个部分的随机成分x＝如果有n个观测值,则样本数据集记为X,即

假定成分数据集X有近似零值,则对应的探测范围记为E＝[e ij]n×D,其中e ij是x ij的探测范围。记R j⊂{1,2,…,n}为第j(j∈{1,2,…,D})个成分中有近似零值的行指标,则O j＝{1,2,…,n}\R j表示第j个成分的剩余行指标,即

首先,我们通过乘法替换方法[56]对近似零值进行初始化,近似零值替换值为探测范围的65%,则X表示替换后的数据集。记随机成分x的ilr坐标为z＝ilr(x)＝[z 1,z 2,…,z D－1]＝[z 1,z－1],其中z－1表示z中除了第一个成分外的剩余成分。因此初始化数据集X变换为实数数据集Z＝[z ij]n×(D－1),其中Z的每行是X中对应成分数据的ilr坐标。对于探测范围E中的元素e i 1,近似零值x i 1＜e i 1通过ilr变换则会导致缺失数据z i 1,该值小于ψi 1,其中

在本节提出的方法中,未知数据z i 1(i∈R 1)通过如下条件期望插补

其中z i,－1表示Z中除第一列外的矩阵中的第i行,条件密度函数f(z 1|z－1＝z i,－1)可以通过如下公式计算

无论多元随机变量z的分布如何,密度函数f(z 1,z－1＝z i,－1)都可以通过多元高斯核密度估计[88]。在这里我们假定不同的坐标方向具有相同的窗宽h,则

窗宽h通过公式h＝σ给出[82],其中σ2＝tr代表矩阵Var(z)的迹。

基于公式(3.1.2)和(3.1.3)可以得到

通过公式(3.1.4)中的条件密度函数,公式(3.1.1)可以表示为

因为

和

其中φ(·)和Φ(·)分别是标准正态分布的密度和分布函数,从而公式(3.1.5)可以化简为

因此,未知数据z i 1可以通过公式(3.1.6)插补。对于ilr坐标,因为d(z i,－1,z k,－1)＝d a(x i,－1,x k,－1),则z i 1的插补值取决于子成分x i,－1与x k,－1之间的Aitchison距离,其中x i,－1与x k,－1分别表示x i与x k中除了第一个成分外的剩余成分。

性质3.1.1　公式(3.1.6)中E(z 1|z－1＝z i,－1,z 1＜ψi 1)插补值有如下性质:

(1)插补值低于探测范围,即E(z 1|z－1＝z i,－1,z 1＜ψi 1)＜ψi 1。

(2)当x中除了第一个成分外的剩余成分任意置换时,插补值不变。

(3)在正交基{e 2,e 3,…,e D－1}改变的情况下插补值是不变的。性质3.1.1是显然的。根据z＝xΨT可以得到

ΨTΨ)＝tr(Var(x)G D),其中G D＝I D－J D,I D是单位矩阵,J D是元素全为1的矩阵[22]。因此所有潜在的元素d(z i,－1,z k,－1),h,z k 1和ψi 1在置换和正交基改变的情况下是不变的,从而公式(3.1.6)中插补值是不变的。

性质3.1.1中(2)和(3)指出E(z 1|z－1＝z i,－1,z 1＜ψi 1)满足置换不变性和正交基改变不变性,但是E(z l|z－l＝z i,－l,z l＜ψil)(l＝2,…,D－1)有可能不满足这两个性质,例如,当x中除了第l个成分的剩余成分任意置换时,z kl有可能改变。为了替换x中第l个成分的近似零值,定义置换后的成分数据为x(l)＝＝x l,x 1,…,x l－1,(x l＋1,…,x D),则ilr坐标为z(l)＝ilr(x(l))＝对应的ilr数据集记为Z(l)＝根据公式(3.1.6),成分数据集X的第l个成分的第i行的近似零值产生的未知数据可以通过如下公式插补。(www.daowen.com)

其中

类似于基于ilr坐标的修正EM算法[59],本节所提出的方法的具体步骤如下:

Step 1:按照成分数据集中每个成分的近似零值比率从大到小进行排序,进而对成分进行置换。为了减小误差,有更多零值的成分应该被放在第一列。不失一般性,假定所有成分已经排序,|R 1|≥|R 2|≥…≥|R D|,其中表示R j(j＝1,2,…,D)的元素个数。

Step 2:通过乘法替换方法对近似零值进行初始替换。

Step 3:设置l＝1。

Step 4:基于公式(3.1.7)替换未知数据

Step 5:基于公式(2.2.6)对更新后的数据集的每行做ilr逆变换。

Step 6:对每个l＝2,3,…|C|执行Steps 4～5,其中C＝{j:j∈{1,2,…,D},|R j|≠0}是包含至少一个近似零值的成分指标。

Step 7:重复Steps 3～6直到上一步与下一步迭代的成分数据集方差矩阵之间的欧氏距离小于一个确定的临界值。

Step 8:对替换后的成分数据集的成分按原始顺序进行排序。

如果数据集X＝闭合为常数,则替换后数据集＝通过以上算法获得。否则,我们应该对替换后值通过如下表达式进行尺度化[62]

其中是尺度化后值,x ik是成分数据集X中第i行第k列原始观测到的元素,是中对应的替换值。

为了验证本节所提出的方法的有效性,我们通过模拟分析和实例分析来比较该方法(mult K)与乘法替换法(mult R)、Kaplan-Meier法、乘法对数正态替换方法(mult LN)、基于alr坐标的修正EM算法(alr EM)、基于ilr坐标的稳健修正EM算法(ilr EM)。给定原始成分数据集X,该数据集中没有近似零值。我们设置低于探测范围以下的值为零,通过替换方法替换后的成分数据集记为X*。考虑两种评价指标:标准残差平方和(STRESS)[56]与方差矩阵之间的相对差异(RDVM)[60]。对于任意随机成分x＝x 1,x 2,…,x D

(),方差矩阵[19]定义为

其中方差矩阵中元素是含有D个部分成分x中任意两个部分i与j之间的对数比率方差。

记原始数据集X和插补后数据集X*的方差矩阵分别为T＝[t ij]D×D和T*＝两种评价指标STRESS和RDVM分别定义为

其中x i是数据集X的第i行。两种评价指标STRESS和RDVM分别表示距离差异和方差差异。