基于多元玻璃物证的非平衡似然比比对检验模型建立及评估

Establishment and Evaluation of Non-Equilibrium Likelihood Ratio Comparison Test Model Based on Multi-Element Glass Physical Evidence

摘要: 对在犯罪现场和嫌疑人身上发现的证据特征进行测量和评估是刑事调查中的重要部分。Aitken等在2006年提出了样本平衡的多元正态分布三层次证据似然比模型。基于玻璃数据，本文考虑了一种非平衡类型的多元正态三层次证据：包含7种元素的玻璃碎片数据58组，其中每组包含的玻璃样本数量不同(样本非平衡性)；共有396个玻璃样本，每个样本都被测量4次。为降低数据维度，利用图形模型(graphical models)找出具有高度相关的变量。假定数据服从多元正态分布，建立样本非平衡的多元正态分布三层次证据似然比模型。在58对同一来源对比中，$LR>1$ 的有20对，假阴性率约为34.48%；在1653对不同来源对比中，$LR>1$ 的有78对，假阳性率约为4.7%。对于司法机关来说，评价的假阳性率更为重要，因为假阳性结果是对应于对无辜者的证据。本文建立的基于图形模型的非平衡似然比建模为司法科学中解决非平衡数据提供了一种新的方案。

Abstract: In criminal investigation, the measurement and evaluation of the characteristics of the evidence found at the crime scene and on the suspect is an important part. Aitken et al. proposed a multivariate normal distribution three-level evidence likelihood ratio model with sample balance in 2006. Based on the glass data, we consider a non-equilibrium type of multivariate normal three-level evidence: 58 sets of glass fragment data containing 7 elements, each of which contains a different number of glass samples (sample non-equilibrium). There were 396 glass samples, each of which was measured four times. To reduce the dimensionality of the data, graphical models are used to find variables with high correlation. In the 58 pairs from the same source, there were 20 pairs $LR>1$ , and the false negative rate was about 34.48%. In the comparison of 1653 pairs from different sources, there were 78 pairs $LR>1$ , and the false positive rate was about 4.7%. For the prosecution, the false positive rate of the evaluation is more important because false positive results correspond to evidence against an innocent person. The non-equilibrium likelihood ratio modeling based on graphical model constructed in this paper provides a new scheme for solving non-equilibrium data in forensic science.

1. 引言

基于贝叶斯推理框架的似然比法[1]-[3] (likelihood ratio, LR)是目前世界公认的法庭科学领域证据价值评价方法之一，它基于先验知识和证据信息对某一事件发生的可能性进行逻辑推理和概率估计，自提出至今已被广泛应用在多个物证鉴定专业领域。例如玻璃碎片分析[4]、枪弹痕迹[5]、墨水笔迹[6]、毒品检验[7]等。

与国外相比，目前国内似然比的研究和应用仍处于起步阶段。虽然已有对某类微量物证进行似然比模型建立和计算的研究[8]，但只是针对单变量的数据。对于更为复杂的、包含多个变量的微量物证，目前还没有相关的研究。

本文以玻璃物证为分析对象，基于多种元素玻璃检验数据，提出了一种评估多变量证据的似然比方法，该方法是对Aitken等[9]提出的样本平衡三层次正态分布证据似然比模型的扩展。为了避免“维度诅咒”，利用图形模型(Graphical model)找出具有高度相关的元素，使得高维数据集变为相互独立的几个低维数据集的乘积，从而尽可能保留原始数据集中的有用信息。

利用两类错误率来评估样品比对检验的结果是法庭科学领域常用的方法[10]。两类错误率是指假阴性率(False negative rate)和假阳性率(False positive rate) [11]来评价LR模型的判断效果。假阴性率是指所有来自同一来源的样本LR结果中，出现$\text{LR}<1$ 的样本个数占比；假阳性率则是指所有来自不同来源的样本LR结果中，出现$\text{LR}>1$ 的样本个数占比。通过两类错误率可以直观地对LR模型的好坏作一个评价。

2. 测量方法

玻璃在日常用品中被广泛使用，使得玻璃成为犯罪现场经常遇到的重要法医痕迹证据。在玻璃的生产过程中，许多不同的元素被纳入熔融混合物中，这些主要成分是硅、钠、钙、镁和钾等的氧化物。采用SEM-EDX法[12] [13]对玻璃样品的元素进行测定。先用树脂包埋法处理需要测量的玻璃，然后利用专业的扫描电子显微镜(SEM-EDX设备)分析玻璃样品的SEM图像和获得的玻璃样品的光谱，最终测定了七种($p=7$ )氧化物的含量，分别是氧化钠(NaO)、氧化镁(MgO)、氧化铝(AlO)、氧化硅(SiO)、氧化钾(KO)、氧化钙(CaO)和氧化铁(FeO)。

2.1. 基础数据排版规范的完整性

考虑$s=58$ 组玻璃样本，每组样本包含不同的样本数量$n_i\left(i=1,\cdots ,s\right)$ (样本非平衡)，$M={\displaystyle \sum }_{i=1}^s{n}_i=396$ ，每个样本独立重复测量$m=4$ 次。这样分组–组内样本–测量次数就构成了三层次框架。为使数据更符合多元正态分布，对数据进行对数变换(${\text{log}}_{10}\left(L\right)$ ，$L$ 为测定的元素含量)。以上所有已被测定的数据从《Statistical Analysis in Forensic Science》获得。

2.2. 样本数据和检材数据

假设在犯罪现场收集到的$G_1$ 个样本证据为样本数据，这些数据的来源假定是已知的，对它们分别进行$N_1$ 次独立重复测量；假设在嫌疑人身上收集到的$G_2$ 个样本证据为检材数据，这些数据的来源假定是未知的，对它们分别进行了$N_2$ 次的独立重复测量。

3. 建立物证比对似然比模型

Aitken等在2006年提出的多元正态分布三层次证据似然比模型，考虑的是样本平衡的数据，即每组都包含相同数量的样本。但在实际法医工作中，更多地会出现样本非平衡的情况，即收集到的每组数据，包含不同数量的样本。如果依然采用样本平衡下的似然比模型，需要对每组数据进行部分删减，毫无疑问会丢失大量有用的信息，以及对后续的似然比结果的准确性产生影响。考虑到这种情况，本文将提出样本非平衡的多元正态分布三层次证据似然比模型。相较于Aitken的模型，该模型对正态分布假设和参数估计部分进行了改进，使其体现出样本非平衡这一特性，从而更好地保留了数据中的信息。

3.1. 正态分布假设

现从基础数据中随机抽样，记第$i$ 组的第$j$ 个样本第$k$ 次测量的证据向量为${\text{X}}_{ijk}={\left(x_{ijk1},\cdots ,x_{ijkp}\right)}^{\text{T}};$ $i=1,\cdots ,s;j=1,\cdots ,n_i;k=1,\cdots ,m$ 。假设三层次数据服从正态分布，具体情形如下：

λ 测量层次：由${}_{\overline{\text{X}}}=\left(1/m\right){\displaystyle \sum }_{k=1}^m{\text{X}}_{ijk}$来估计第$i$ 组的第$j$ 个样本的测量均值向量${\theta }_{ij}$ ，假设所有样本的重复测量的协方差矩阵为$U$ ，则${\text{X}}_{ijk}$ 的分布为$\left({\text{X}}_{ijk}|{\theta }_{ij},U\right)~N\left({\theta }_{ij}|U\right);i=1,\cdots ,s;j=1,\cdots ,n_i;k=1,\cdots ,m$ ；

λ 组内层次：由${}_{\overline{\text{X}}}=\left(1/n_i\right){\displaystyle \sum }_{j=1}^{n_i}{}_{\overline{\text{X}}}$来估计第i组的组内均值向量${\sigma }_i$ ，假设组内协方差矩阵为$V$ ，则测量均值向量${\theta }_{ij}$ 的分布为$\left({\theta }_{ij}|{\sigma }_i,V\right)~N\left({\sigma }_i|V\right);i=1,\cdots ,s;j=1,\cdots ,n_i$ ；

λ 组间层次：由${}_{\overline{\text{X}}}=\left(n_i/M\right){\displaystyle \sum }_{i=1}^s{}_{\overline{\text{X}}}$来估计总体均值向量$\mu$ ，假设组间协方差矩阵为$W$ ，则组内均值向量${\sigma }_i$ 的分布为$\left({\sigma }_i|\mu ,W\right)~N\left(\mu |V\right);i=1,\cdots ,s$ 。

假设检材数据和样本数据分别记为$y_1,y_2$ 。为方便以下表述，引入指标$d$ ：$d=1$ 代表检材数据，$d=2$ 代表样本数据。记指标$d$ 的第$j$ 个样本第$k$ 次测量的证据向量为$y_{djk}={\left(y_{djk1},\cdots ,y_{djkp}\right)}^{\text{T}};$ $j=1,\cdots ,G_d;k=1,\cdots ,N_d$ 。此时，记指标$d$ 的第$j$ 个样本的重复测量的样本均值向量为${\overline{y}}_{dj}=\left(1/N_d\right){\displaystyle \sum }_{k=1}^{N_d}{y}_{djk}$ 。记指标$d$ 的重复测量的样本均值向量为${\overline{y}}_d=\left(1/G_d\right){\displaystyle \sum }_{j=1}^{G_d}{\overline{y}}_{dj}$ 。现定义加权均值$\overline{y}$ ，则$\overline{y}=\left(G_1{N}_1{\overline{y}}_1+G_2{N}_2{\overline{y}}_2\right)/t$ ，其中$t=G_1{N}_1+G_2{N}_2$ ，加权均值$\overline{y}$ 的分布为：$\left(\overline{y}|\mu ,U,V,W\right)~N\left(\mu ,{\left(G_d{N}_d\right)}^{-1}U+N_d^{-1}V+W\right)$ 。

3.2. 参数估计

按照经典的极大似然理论[14]，总体均值向量$\mu$ 由${}_{\overline{\text{X}}}$来估计。

利用$\left\{{\text{X}}_{ijk}\right\}$ 估计测量层次的协方差矩阵$U$ ，其表达式如下：

$\widehat{U}=\frac{S_U}{M\left(m-1\right)}$ (1)

其中，$S_U={\displaystyle \sum }_{i=1}^s{\displaystyle \sum }_{j=1}^{n_i}{\displaystyle \sum }_k^m\left({\text{X}}_{ijk}-{}_{\overline{\text{X}}}\right){\left({\text{X}}_{ijk}-{}_{\overline{\text{X}}}\right)}^{\text{T}}$。

利用$\left\{{\text{X}}_{ijk}\right\}$ 估计组内层次的协方差矩阵$V$ ，其表达式如下：

$\widehat{V}=\frac{S_V}{M-s}-\frac{\widehat{U}}{m}$ (2)

其中，$S_V={\displaystyle \sum }_{i=1}^s{\displaystyle \sum }_{j=1}^{n_i}\left({}_{\overline{\text{X}}}-{}_{\overline{\text{X}}}\right){\left({}_{\overline{\text{X}}}-{}_{\overline{\text{X}}}\right)}^{\text{T}}$。

利用$\left\{{\text{X}}_{ijk}\right\}$ 估计组间层次的协方差矩阵$W$ ，其表达式如下：

$\widehat{W}=\frac{S_W}{s-1}-\frac{\widehat{U}}{m}-\frac{\widehat{V}}{m}$ (3)

其中，$S_w={\displaystyle \sum }_{i=1}^s{n}_i\left({}_{\overline{\text{X}}}-{}_{\overline{\text{X}}}\right){\left({}_{\overline{\text{X}}}-{}_{\overline{\text{X}}}\right)}^{\text{T}}$。

3.3. 非平衡多元正态性假设的似然比模型

从犯罪现场采集的证据样本为$p_1$ ，从嫌疑人身上采集的证据样本为$p_2$ ，$E$ 为由$p_1$ 和$p_2$ 组成的证据样本总体。记控方假设(Prosecution hypothesis)和辩方假设(Defense hypothesis)分别为H_p、H_d [15]，具体意义如下：

H_p：$P_1=P_2$ ，二者是相同来源，认为嫌疑人和犯罪现场存在相关。

H_d：$P_1\ne P_2$ ，二者不是相同来源，认为嫌疑人和犯罪现场不存在相关。

对这两种竞争假设条件下获得证据的概率比进行计算，即似然比为$\text{LR}=P\left(E|H_p\right)/P\left(E|H_d\right)$ 。

基于控方假设H_p和辩方假设H_d，考虑连续性三层次非平衡的正态分布模型样品，参考文献[9]，可得LR值的计算公式如下：

$\text{LR}=\frac{f\left(E|H_p\right)}{f\left(E|H_d\right)}=\frac{f\left(y_1,y_2|\mu ,U,V,W\right)}{f\left(y_1|\mu ,U,V,W\right)f\left(y_2|\mu ,U,V,W\right)}$ (4)

其中，f表示多元正态分布的概率密度函数。分子计算如下：

$\begin{array}{l}f\left(y_1,y_2|\mu ,U,V,W\right)\\ ={\left|2\pi \left(V+W\right)\right|}^{-1/2}{\left|2\pi U\right|}^{-T/2}\times {\left|2\pi {\left\{t{U}^{-1}+{\left(V+W\right)}^{-1}\right\}}^{-1}\right|}^{1/2}\times \text{exp}\left\{-\frac{1}{2}\left(H_1+H_2\right)\right\}\\ \times \text{exp}\left\{-\frac{1}{2}{\displaystyle \sum }_{j=1}^{G_1}\text{tr}\left(S_{1j}{U}^{-1}\right)-\frac{1}{2}{\displaystyle \sum }_{j=1}^{G_2}\text{tr}\left(S_{2j}{U}^{-1}\right)-\frac{1}{2}{N}_1\text{tr}\left(S_1{U}^{-1}\right)-\frac{1}{2}{N}_2\text{tr}\left(S_2{U}^{-1}\right)\right\}\end{array}$

其中：

$D_1=N_1{G}_1^{-1}U;D_2=N_2{G}_2^{-1}U;H_1=\left({\overline{y}}_1-{\overline{y}}_2\right){\left(D_1+D_2\right)}^{-1}{\left({\overline{y}}_1-{\overline{y}}_2\right)}^{\text{T}}$

$H_2=\left(\overline{y}-\mu \right){\left(T^{-1}U+\left(V+W\right)\right)}^{-1}{\left(\overline{y}-\mu \right)}^{\text{T}};S_d={\displaystyle \sum }_{j=1}^{N_d}\left({\overline{y}}_{dj}-{\overline{y}}_d\right){\left({\overline{y}}_{dj}-{\overline{y}}_d\right)}^{\text{T}};$

$S_{dj}={\displaystyle \sum }_{j=1}^{N_d}\left(y_{djk}-{\overline{y}}_{dj}\right){\left(y_{djk}-{\overline{y}}_{dj}\right)}^{\text{T}};d=1,2$

分母计算如下：

$\begin{array}{l}f\left(y_1|\mu ,U,V,W\right)\\ ={\left|2\pi U\right|}^{-\left(N_1{G}_1\right)/2}{\left|2\pi \left(V+W\right)\right|}^{-1/2}{\left|2\pi \left[D_1^{-1}+{\left(V+W\right)}^{-1}\right]\right|}^{-1/2}\\ \times \text{exp}\left\{-\frac{1}{2}{\displaystyle \sum }_{j=1}^{G_1}\text{tr}\left(S_{1j}{U}^{-1}\right)-\frac{1}{2}{N}_1\text{tr}\left(S_1{U}^{-1}\right)\right\}\text{exp}\left\{-\frac{1}{2}\left({\overline{y}}_1-\mu \right){\left(D_1+V+W\right)}^{-1}{\left({\overline{y}}_1-\mu \right)}^T\right\}\end{array}$

同理，可得$f\left(y_2|\mu ,U,V,W\right)$ 的表达式：

$\begin{array}{l}f\left(y_2|\mu ,U,V,W\right)\\ ={\left|2\pi U\right|}^{-\left(N_2{G}_2\right)/2}{\left|2\pi \left(V+W\right)\right|}^{-1/2}{\left|2\pi \left[D_2^{-1}+{\left(V+W\right)}^{-1}\right]\right|}^{-1/2}\\ \times \text{exp}\left\{-\frac{1}{2}{\displaystyle \sum }_{j=1}^{G_2}\mathrm{tr}\left(S_{2j}{U}^{-1}\right)-\frac{1}{2}{N}_2\mathrm{tr}\left(S_2{U}^{-1}\right)\right\}\text{exp}\left\{-\frac{1}{2}\left({\overline{y}}_2-\mu \right){\left(D_1+V+W\right)}^{-1}{\left({\overline{y}}_2-\mu \right)}^T\right\}\end{array}$

将两者相乘，得到(4)式的分母。

将以上得到的结果相应地代入(4)式的分子分母，经过简化运算，可以得到LR值的最简表达式如下：

$\text{LR}=\frac{ {\left|V+W\right|}^{1/2}{\left|t{U}^{-1}+{\left(V+W\right)}^{-1}\right|}^{-1/2}\exp \left\{-\frac{1}{2}\left(H_1+H_2\right)\right\}}{{\left|D_1^{-1}+{\left(V+W\right)}^{-1}\right|}^{-1/2}{\left|D_2^{-1}+{\left(V+W\right)}^{-1}\right|}^{-1/2}\text{exp}\left\{-\frac{1}{2}\left(H_3+H_4\right)\right\}}$ (5)

4. 图形模型

图形模型是一种基于偏相关矩阵和图论的降维技术。该模型具有很好的变量可解释性，同时减少了需要估计的参数数量。在图中，每个变量代表一个节点，具有高度相关的变量由一条线(边)连接。几个具有相关性的变量联结在一起称为团。而不同变量间的相关性由偏相关矩阵来确定。在本文中，偏相关矩阵由组间协方差矩阵和组内协方差矩阵之和$\left(V+W\right)$ 确定。通过先对$\left(V+W\right)$ 求逆，再进行缩放，可得到偏相关矩阵[16]。一般地，偏相关矩阵是正定的；如果仅是半正定，可通过重新参数化变成正定[17]。

在实际数据分析时，需要设定一个相关系数的阈值$R_0$ ，$0<R_0<1$ 。对于$\Omega$ 中的相关系数$R_{ab}\left(a,b=1,\cdots ,p,a\ne b\right)$ ，若$\left|R_{ab}\right|\ge R_0$ ，则保留该相关系数；否则令$R_{ab}=0$ 。$R_0$ 取值不同，得到的偏相关矩阵(记为${\Omega }_{R_0}$ )也不同。根据不同的偏相关矩阵会产生不同的图形模型。

利用图形模型，在计算似然比时，可以把高维的数据集，变为几个低维的数据集的乘积。计算公式如下：

$\mathrm{LR}\left(C_i|S_i\right)=\frac{\prod \mathrm{LR}\left(C_i\right)}{\prod \mathrm{LR}\left(S_i\right)}$ (6)

其中，$C_i$ 为图形模型中的第$i$ 个团，即所有节点相互连接的变量子集。第$i$ 个团的所有分隔符的集合用$S_i$ 表示，即两个团相交处的一个节点或一组节点。

5. 实验与分析

5.1. 实验流程

在第2.1节，给出了多种化合物玻璃数据的具体信息；假设其中被测量的7个变量在经过对数变换后分布为7维的多元正态分布，7维向量的分量为log₁₀(NaO)、log₁₀(MgO)、log₁₀(AlO)、log₁₀(SiO)、log₁₀(KO)、log₁₀(CaO)、log₁₀(FeO)。

5.2. 参数结果

基于基础数据，利用第3.1节的公式依次计算，得到参数$\mu$ 的估计：

$\widehat{\mu }={\left(-0.7095,-2.0796,-2.2549,-0.1408,-2.8842,-1.2408,-4.6718\right)}^{\text{T}}$

参数U的估计：

$\widehat{U}=\left(\begin{array}{ccccccc}2.029\times {10}^{-4}& 3.051\times {10}^{-4}& -2.177\times {10}^{-4}& 1.813\times {10}^{-5}& -1.923\times {10}^{-4}& \cdots 2.799\times {10}^{-4}& 2.980\times {10}^{-4}\\ 3.051\times {10}^{-4}& 1.711\times {10}^{-1}& -1.387\times {10}^{-3}& 3.387\times {10}^{-4}& -1.633\times {10}^{-3}& -3.671\times {10}^{-3}& -3.897\times {10}^{-3}\\ -2.177\times {10}^{-4}& 1.387\times {10}^{-3}& -1.528\times {10}^{-1}& 2.984\times {10}^{-4}& 1.810\times {10}^{-3}& -6.095\times {10}^{-6}& 1.848\times {10}^{-3}\\ 1.813\times {10}^{-5}& 3.387\times {10}^{-4}& -2.984\times {10}^{-4}& 8.009\times {10}^{-4}& 2.275\times {10}^{-3}& 1.123\times {10}^{-3}& 6.104\times {10}^{-4}\\ -1.923\times {10}^{-4}& -1.633\times {10}^{-3}& -1.810\times {10}^{-3}& 2.275\times {10}^{-3}& 2.531\times {10}^{-1}& 1.203\times {10}^{-3}& -2.318\times {10}^{-2}\\ 2.799\times {10}^{-4}& -3.671\times {10}^{-3}& -6.059\times {10}^{-6}& 1.123\times {10}^{-3}& 1.203\times {10}^{-3}& -3.352\times {10}^{-2}& -4.165\times {10}^{-4}\\ 2.980\times {10}^{-4}& -3.897\times {10}^{-3}& 1.848\times {10}^{-6}& 6.104\times {10}^{-4}& -2.318\times {10}^{-2}& 4.165\times {10}^{-3}& 3.104\times {10}^{-4}\end{array}\right)$

参数V的估计：

$\widehat{V}=\left(\begin{array}{ccccccc}3.140\times {10}^{-3}& 1.966\times {10}^{-2}& -8.823\times {10}^{-4}& 7.461\times {10}^{-6}& -1.109\times {10}^{-2}& 1.989\times {10}^{-2}& 2.980\times {10}^{-3}\\ 1.966\times {10}^{-2}& 6.562\times {10}^{-1}& -1.030\times {10}^{-1}& -7.660\times {10}^{-4}& -8.508\times {10}^{-2}& 2.584\times {10}^{-1}& -3.897\times {10}^{-2}\\ 8.823\times {10}^{-4}& -1.030\times {10}^{-1}& 3.317\times {10}^{-1}& -1.609\times {10}^{-3}& 1.016\times {10}^{-1}& -1.965\times {10}^{-2}& 1.848\times {10}^{-1}\\ 7.461\times {10}^{-6}& -7.660\times {10}^{-4}& -1.609\times {10}^{-3}& 1.323\times {10}^{-3}& 6.330\times {10}^{-3}& 1.656\times {10}^{-3}& 6.104\times {10}^{-3}\\ -1.109\times {10}^{-2}& -8.508\times {10}^{-2}& 1.016\times {10}^{-1}& 6.330\times {10}^{-3}& 8.868\times {10}^{-1}& -7.847\times {10}^{-2}& -2.318\times {10}^{-2}\\ 1.989\times {10}^{-2}& 2.584\times {10}^{-1}& -1.965\times {10}^{-2}& 1.656\times {10}^{-3}& -7.847\times {10}^{-2}& 3.148\times {10}^{-1}& 4.165\times {10}^{-2}\\ 5.381\times {10}^{-3}& 3.062\times {10}^{-2}& -1.221\times {10}^{-1}& 3.363\times {10}^{-3}& 7.784\times {10}^{-2}& 2.919\times {10}^{-2}& 1.1110\end{array}\right)$

参数W的估计：

$\widehat{W}=\left(\begin{array}{ccccccc}1.172& 3.548& 3.234& 0.229& 3.971& 2.559& 7.767\\ 3.548& 20.26& 7.913& 0.665& 5.661& 11.06& 26.29\\ 3.234& 7.913& 13.89& 0.640& 16.77& 5.132& 19.42\\ 0.229& 0.665& 0.640& 0.049& 0.812& 0.444& 1.516\\ 3.971& 5.661& 16.77& 0.812& 27.46& 4.056& 24.69\\ 2.559& 11.06& 5.132& 0.444& 4.056& 11.37& 17.03\\ 7.767& 26.29& 19.42& 1.516& 24.69& 17.03& 58.4\end{array}\right)$

5.3. 建立图形模型

根据图形模型方法，先求协方差矩阵估计$\left(V+W\right)$ 的逆，再缩放得到逆缩放矩阵$\Omega$ ，如表1所示。

Table 1. Scaled inverse matrix $\Omega$

表1. 缩放逆矩阵$\Omega$

	logNaO	logMgO	logAlO	logSiO	logKO	logCaO	logFeO
logNaO	1.000	0.0019	−0.4119	−0.6394	−0.0413	−0.5071	−0.5820
logMgO		1.000	−0.2569	−0.0347	0.4348	−0.3062	−0.3720
logAlO			1.000	0.0956	−0.6360	0.1500	0.3807
logSiO				1.000	−0.0740	0.2653	0.0172
logKO					1.000	0.1529	−0.2288
logCaO						1.000	0.1882
logFeO							1.000

取阈值$R_0=0.5$ ，对于$\Omega$ 来说，$\left|R\right|<0.5$ 的值全部设定为0，可以得到第一个偏相关矩阵${\Omega }_{0.5}$ ：

${\Omega }_{0.5}=\left(\begin{array}{ccccccc}1.0000& \text{*}& \text{*}& -0.6394& \text{*}& -0.5071& -0.5820\\ & 1.0000& \text{*}& \text{*}& \text{*}& \text{*}& \text{*}\\ & & 1.0000& \text{*}& -0.6360& \text{*}& \text{*}\\ & & & 1.0000& \text{*}& \text{*}& \text{*}\\ & & & & 1.0000& \text{*}& \text{*}\\ & & & & & 1.0000& \text{*}\\ & & & & & & 1.0000\end{array}\right)$

将偏相关矩阵${\Omega }_{0.5}$ 中的相关系数绝对值从大到小排序，可以确定以下的图形模型(见图1)。

Figure 1. Graphical model of partial correlation matrix ${\Omega }_{0.5}$

图1. 偏相关矩阵${\Omega }_{0.5}$ 的图形模型

由图1得到，CaO、NaO、SiO、FeO具有相关性，AlO、KO具有相关性，化合物MgO与其他元素没有相关性。进一步，由图1可以确定基于偏相关矩阵${\Omega }_{0.5}$ 的LR值计算公式：

$\mathrm{LR}=\frac{\mathrm{LR}\left(\text{Ca},\text{Si}\right)\mathrm{LR}\left(\text{Si},\text{Na}\right)\mathrm{LR}\left(\text{Fe},\text{Na}\right)\mathrm{LR}\left(\text{Al},\text{K}\right)\mathrm{LR}\left(\text{Mg}\right)}{\mathrm{LR}\left(\text{Na}\right)\mathrm{LR}\left(\text{Na}\right)}$ (7)

若取阈值$R_0=0.4$ ，对于$\Omega$ 来说，$\left|R\right|<0.4$ 的值全部设定为0，可以得到第二个偏相关矩阵${\Omega }_{0.4}$ ：

${\Omega }_{0.4}=\left(\begin{array}{ccccccc}1.0000& \text{*}& -0.4119& -0.6394& \text{*}& -0.5071& -0.5820\\ & 1.0000& \text{*}& \text{*}& 0.4348& \text{*}& \text{*}\\ & & 1.0000& \text{*}& -0.6360& \text{*}& \text{*}\\ & & & 1.0000& \text{*}& \text{*}& \text{*}\\ & & & & 1.0000& \text{*}& \text{*}\\ & & & & & 1.0000& \text{*}\\ & & & & & & 1.0000\end{array}\right)$

将偏相关矩阵${\Omega }_{0.4}$ 中的相关系数绝对值从大到小排序，可以确定以下的图形模型(见图2)。

Figure 2. Graphical model of partial correlation matrix ${\Omega }_{0.4}$

图2. 偏相关矩阵${\Omega }_{0.4}$ 的图形模型

由图2得到，7种元素依次具有相关性。再由图2，可以确定基于偏相关矩阵${\Omega }_{0.4}$ 的LR值计算公式：

$\mathrm{LR}=\frac{\mathrm{LR}\left(\text{Ca},\text{Na}\right)\mathrm{LR}\left(\text{Si},\text{Na}\right)\mathrm{LR}\left(\text{Fe},\text{Na}\right)\mathrm{LR}\left(\text{Na},\text{Al}\right)\mathrm{LR}\left(\text{Al},\text{K}\right)\mathrm{LR}\left(\text{K},\text{Mg}\right)}{\mathrm{LR}\left(\text{Na}\right)\mathrm{LR}\left(\text{Na}\right)\mathrm{LR}\left(\text{Na}\right)\mathrm{LR}\left(\text{Al}\right)\mathrm{LR}\left(\text{K}\right)}$ (8)

5.4. 同一来源的玻璃数据比对

为进行比对实验，现在从基础数据中抽取样本数据和检材数据，将剩余的数据用于计算参数。具体操作如下，从$s=58$ 组中随机抽取一组，选取该组的任意一个样本作为检材数据$\left(N_2=1\right)$ ，该组余下的样本作为样本数据。本文做实验时，仅设定选取第一个样本作为检材数据；剔除已被选定的一组，将余下的57组数据用作参数估计。先利用第5.3节给出的公式(7)或公式(8)确定LR值计算公式，再进行比对检验。最终每个图形模型都输出了58个比对检验结果，根据不同数量级，输出LR值计数统计结果见表2。

Table 2. Comparison of two graphical models when recovered data and control data from the same source

表2. 样本数据和检材数据同源时两种图形模型的比较

似然比	$\left|R\right|<0.5$	$\left|R\right|<0.4$
<1	20	20
1~10	1	1
10~10^2	1	1
10^2~10^3	6	6
10^3~10^4	15	15
10^4~10^5	15	15
>10^5	0	0
所有样本	58	58
真阳性(%)	65.52	65.52
假阴性(%)	34.48	34.48

5.5. 不同来源的玻璃数据比对

具体操作如下，从58组数据中随机抽取两组，从抽取的第一组中，随机地抽取一个样本作为检材数据$\left(N_2=1\right)$ ，将第二组全体作为样本数据。剔除已被抽取的两组，将余下的56组数据用作参数估计。随机抽取两组的组合数为${\mathrm{C}}_{58}^2=1653$ ，最终输出了1653个比对检验结果，计数统计结果见表3。

Table 3. Comparison of two graphical models when recovered data and control data from the different source

表3. 样本数据和检材数据不同来源时两种图形模型的比较

似然比	$\left|R\right|<0.5$	$\left|R\right|<0.4$
>1	78	76
1~10^−^1	17	16
10^−^1~10^−^2	18	16
10^−^2~10^−^3	16	20
10^−^3~10^−^4	12	7
10^−^4~10^−^5	17	20
<10^5	1495	1498
所有样本	1653	1653
真阴性(%)	95.28	95.40
假阳性(%)	4.718	4.597

6. 结论

对于同一来源的玻璃碎片，假阴性的比例比较大，有接近35%的样品似然比小于1，并被判别为来自不同来源。对于不同来源的玻璃碎片，假阳性占比非常低，只有4.7%左右的碎片似然比大于1，并被判别为来自同一来源。虽然相较于Aitken和Lucy的方法，两类错误率下降较为明显(在Aitken的方法中，其假阴性率大于40%，假阳性率在7%左右)，但可以发现，假阴性的比例依然较高。

在文献[9]中，对两类错误率的高低情况作出了具体的解释。司法机关会更看重假阳性结果，因为假阳性结果对应于无辜者的判别结果。而假阴性对于后续的司法工作影响较小。

在进行数据降维时，分别使用了由$\left|R\right|<0.5$ 或$\left|R\right|<0.4$ 导出的两组图形模型。由两类错误率来看，两组模型基本没有太大的区别。但是考虑到$\left|R\right|<0.4$ 时导出的图形模型，所有元素都存在相关性。因此，考虑使用$\left|R\right|<0.4$ 时的图形模型。

对于司法机关来说，似然比强调的是物证检验结果本身的概率之比，其数值大小反映的是物证支持两个假设的程度，而不是指来源可能性之比。因此，在司法实践中，要想回答物证是否同源、同源的可能性有多大等关键问题，法官还需要采用似然比与其他证据相互印证、排除有效怀疑的方式等诸多方法[18]进行进一步的评估。

本文描述了一种以非平衡三层次模型的形式评价证据的方法。该方法可适用于高维数据，通过考虑由缩放逆协方差矩阵得到的图形模型降低维度进行似然比计算，最终提供了一个低假阳性率的结果，这在法医科学中是可以接受的。

基金项目

本课题来自公安部物证鉴定中心《基于理化检验数据的似然比证据价值评价方法研究》，是公安部物证鉴定中心基本科研业务费专项资金自主选题项目，项目编号为2022JC12。