章11-2-模式识别方法

§11.5.7．SIMCA方法[13]SIMCA(softindependentmodellingofclassanalogy或者statisticalisolinearmulticategoryanalogy)属于类模型方法，即对每类构造一主成分回归的数学模型，并在此根底上去进行试样的分类。此法在1976年由瑞典学者S.Wold所提出，很快受到普遍的重视，并在化学中得到广泛的应用。如有数据阵其中，n为试样数，m为变量，即维数。对于某一类，主成分回归模型为式中，ai为变量i之均值，A为主成分数，βia是变量i在主成分a上的装载,θak是试样k是关于主成分a的得分,εik是偏差.对于多类，那么主成分回归模型为

其中q表示类。SIMCA方法分类的主要操作步骤可归结为：1.数据标准化假设数据来源不一时，变量间彼此差异可能很大，那么须将数据进行标准化(autoscaling)处理。经过标准化处理的每一变量均值为0，方差为1(见§11.2)。2.主成分数A的测定确定A值的最有效的方法为交叉验证(crossvalidation)法。这种方法的操作过程为：(1)将训练集中的某类分成T组，在分组中要尽可能考虑到试样的代表性。如某类由25个试样所组成，令T=8，即将之分为8组。第一组为1，9，17和25，第二组为2，10，和18，第三组为3，11和19，等，直到第八组，试样为8，16和24。(2)首先将第一组试样从训练集中除去，并设降维的数据矩阵为Y--，试样数为n-.(3)对于Y--，应用前述单类主成分模型公式(11.9)去拟合。拟合中依次令A=0，1，2，…直到M-2或n--2(取决于二者中小者)(4)运用在(3)中建立的数学模型去拟合所除去的试样。此步中，A=0，1，…，并且α,β固定不变。相应于每一A值，计算试样的偏差εik,由此得到这些偏差平方的加和ΔA。(5)将所除去的那组试样重新放回数据阵Y。(vi)由数据阵Y中除去下一组试样，从而得到一新的降维数据阵Y--，回到步(3).假设每一组均被除去一次，那么到第(7)步。(7)对于每一A值，将ΔA加和得到DA，由(DA-1-DA)/n对DA/[n(M-A-1)]作F检验来判断A的重要性，从而确定A值。一旦每类A值确定之后那么可使式(11.10)，即多类主成分模型对于试样的判别能力到达极大。3.主成分模型中β，θ等参数的测定式(11.10)中βia,θak由矩阵Z(q)’Z(q)对角化求得。其中Z(q)为第q类训练集中每一变量减去平均值后所形成的数据阵，Z(q)’为Z(q)的转置矩阵.εik(q)可由Z值减去公式中β和θ的乘积项得到，那么方差可由下式求得：一旦每一类中上述参数求出之后，即可运用主成分回归模型去预测未知试样。4.未知试样测试用主成分回归模型去拟合未知试样Yip.与一般多元回归相同，此时Zi=yip-ai(q)为因变量，βia(q)(a=1,2,…,Aq)为自变量：此处ca为回归系数。表征试样p拟合好坏，可由偏差εip的方差表示：判别试样p是否属于q类，那么用F显著性检验。对于类内试样F检验的公式为：其中，校正因子nq/(nq-Aq-1)的参加是因为当计算该类β和θ时P试样已被包括进去。将F的计算值与临界值[自由度分别为(M-Aq)和(nq-Aq-1)(M-Aq)]相比较，假设F<F临界，那么试样p归入q类，否那么，将试样p拟合于其他类。当试样p拟合于其他类时，即p已不是类内试样，其F显著性检验的公式为：同样将F的计算值与临界值[自由度分别为(M-Aq)和(nq-Aq-1)(M-Aq)]相比较，假设计算值小于临界值，那么将试样p归入其他某q类，否那么，试样p为一新类。5.两类间相似度—非相似度测量用类r中的所有试样去拟合类q主成分模型，那么可得到类间的相似度测量，类r与q间的方差为：同样可以得到S(r,p)2。将式(11.13)所得方差与式(11.11)所得方差相比较，那么可得两类间相似度的测量。6.变量重要性的测量变量在判别中的重要性可由剩余方差(residualvariance)与原始数据的方差相比较而得。假设原始数据经过标准化处理，那么所有变量i的方差相同：由此可得Ui值越大，即剩余方差与原始数据方差的比值愈小，该变量在主成分模型中的作用越小。7.试样相关性测量与变量相类同，即将试样的剩余方差式(11.12)与某一类整个方差式(11.11)相比较(F检验)，其剩余方差越小，该试样与此类的相关性越大。[例]鸢尾花的分类。在30年代，由R.A.Fisher所收集的有关鸢尾花的数据，是后来常被用于模式识别的一套经典的数据。我们取其一局部，即第一类，多刚毛(setose)鸢尾花和第二类，杂色(versicolor)鸢尾花。每类25个试样，两类共50个。对花的描述是：〔1〕花萼长度(x1)；〔2〕花萼宽度(x2)；〔3〕花瓣长度(x3)；〔4〕花瓣宽度(x4)。这50个试样的测试数据及相应的归类示于表11.9。表11.9鸢尾花的测试数据及归类对于原始数据，每个试样按照类分别减去均值。然后，按照类分别进行主成分分析。取一个主成分，那么在主成分回归方程式中βia(装载)和θak(得分)分别示于表11.10和表11.11。表11.10主成分回归方程中的装载(βia)表11.11主成分回归方程中的得分(θak)类1和类2的标准方差分别为S0(1)2=0.0257和S0(2)2=0.0525。试样分别拟合于两类主成份回归模型的标准方差及F检验分别示于表11.12。两类显著性检验F的临界值相同，即均为F(3,69)0.01=4.08。在表中第二和第三列，括弧中数码1或2分别表示用类1模型计算或用类2模型计算所得方差。所以，“最近类〞，即为此试样最有可能所属的类。“第2个最近类〞可依此类推。但，是否属于类1或类2，那么还要通过F检验来确定。如前所述，假假设F的计算值小于F的临界值，那么将试样归于此类，否那么，归于其它类。假设某试样不能归于任一类，那么该试样为新的一类。依此，观察表11.12那么可发现，对于类1模型，试样1-25，只有第23的F计算值〔4.64〕大于4.08〔临界值〕，故试样23不能归于类1；试样26-50可全部正确归于类2。相类同，对于类2模型，试样26-50均可正确地归于类2；而1-25，除了试样23，其它均可正确归于类1。由于试样23的F计算值大于临界值4.08，故亦不能归于类2。试样23为一异常试样，应归于新的一类。在50个试样中，仅有一个试样不能正确归类，所以识别率为98%。表11.12SIMCA方法的分类结果§11.6无管理方法§11.6.1系统聚类分析[14]聚类分析是数理统计中的一种方法，特别适用于试样归属不清楚的情况。它所基于的主要思想是在多维空间中，同类化合物应靠得近些，彼此间的距离小些；相反，不同类的化合物应离得远些，彼此间的距离大些。聚类分析即为如何使相似的试样“聚〞在一起，从而到达分类的目的。聚类分析为无管理方法，其中用得最多的为系统聚类(hierarchicalclusrering)法。系统聚类法的根本思想是首先定义试样之间和类与类之间的距离。在各自成类试样中，将距离最近的两类合并，重新计算新类与其他类间的距离，并按最小距离归类，重复此过程，每次减少一类，直到所有的试样成为一类为止。其聚类过程用图表示，称为聚类图。定义类与类间距离有多种方法，不同的定义就产生了不同的系统聚类分析方法。常用的方法有8种。这些方法是：1.最短距离法此种方法定义Gi与Gj之间的距离为其中dkl是试样xk与xl的距离。也就是说两类之间的距离等于两类中最近试样之间的距离。设某一步将类Gp和Gq合并成Gr，那么类Gi与Gr距离的递推公式为：2最长距离法在此种方法中类与类之间的距离等于两类中最远试样的距离，即显然，最长距离法的递推公式是：3.中间距离法在中间距离法中，类与类之间的距离既不采用两类之间最近的距离，也不采用最远距离，而是采用最远和最近之间的距离。如果在某一步将类Gp和Gq合并为Gr，任一类Gi与Gr的距离的取法可由如右边的三角形说明。如果用最短距离法，Dir等于Dip和Diq中短一边；如果用最长距离法，Di等于Dip和Diq中长的边。而在中间距离法中，Dir既不取Dip，也不取Diq，Dir取的是夹在这两条边之间的三角形中线。由初等几何知道：当β=-1/4时，Dir就是上述的三角形中线。上式即为中间距离法的递推公式。4.重心法该类方法定义两类之间的距离为对应这两类重心之间的距离。对试样来说，每一类的重心即为该类试样的均值。设某一步将Gp与Gq合并成Gr，它们各含有np,nq,nr(nr=np+nq)个试样，那么此种方法的递推公式为：5.类平均法在该种方法中，类Gp与Gq之间的距离表示为：其中，np，nq

为类Gp，Gq中的试样数。就是说类之间的平方距离等于各元素两两之间的平方距离的平均。其递推公式为：6.可变类平均法类平均法的距离递推公式中对于类Gp和Gq之间的距离没有反映进去，可变类平均法将之改进为：其中，β可变，β<17.可变法此种方法的递推公式为：其中，β可变，β<18.方差平方和法这个方法首先由Ward提出，故文献中常称之为Ward法。该方法所基于的思想为方差分析。如某一步中，G1={x1,x2}={1,2},G2={x3,x4}={4.5,6},G3={x5}={8},假设将G1，G2合并，那么两类中所有元素的均值为(1+2+4.5+6)/4=3.375，那么方差平方和为(1-3.375)2+(2-3.375)2+(4.5-3.375)2+(6-3.375)2=15.687。假设将G1与G3合并，那么其方差平方和为(1-3.667)2+(2-3.667)2+(8-3.667)2=28.667。假设将G2和G3合并，那么其方差平方和为6.117。在这三种分类方案中，以G2和G3合并增加的方差平方和为最小，故将G2和G3合并。也就是说，这种方法是将某一步中G1，G2，…Gk类合并成k-1类，而由此增加的方差平方和为最小。该种方法的递推公式为：上述几种方法总的递推公式为：方法不同，式中αp,αq,β,γ的取值也不同。取值方法示于表11.13。其中，ni,np,nq和nr为相应类中试样数。在可变法及可变类平均法中的β可变，分类效果与β取值关系极大。β假设近于1，分类效果不好。通常β取负值。表11.13系统聚类法参数表下面我们以一极为简单的例子来说明聚类分析方法的步骤。如有5个试样，每一试样仅有一个特征：1，2，4.5，6，8。表11.14D0(1)计算试样之间的距离此处选用绝对距离(见§11.4.2)，得表11.14。(2)定义类与类之间的距离如采用最短距离法。(3)逐步归类开始时，5个试样各自成一类，类间距离即为试样间距离，即Dij=dij。选择最小的Dij，此时为类G1和G2，将它们合并成新类G6，并计算与其类间的距离，得表11.15。表11.15D1重复如上过程可得表11.16和11.17，其归类结果表示于图11.9。表11.16D2表11.17D3图11.9系统聚类图最后我们给出一实例：在作者的实验室中曾经测定过血中Sr,Cu,Mg和Zn等4种元素。在前面，此例虽然已被引用过，但是为方便读者，此处仍给出原始测定结果(表11.18)，其中，1—13为冠心病者，14—26为健康人。这4种元素分别与变量x1,x2,x3和x4相对应，即每一样品为由上述变量所构成的4维空间中的一个点。

试样间距离采用塔尼莫特(Tanimoto)系数法；类间距离采用方差平方和(Ward)法，其聚类结果示于图11.10。图11.10血样4元素系统聚类图*ID为试样序号表11.18血中4种元素测定结果(µg/ml)1)测试结果为二次测定平均值。§11.6.2．最小生成树最小生成树法是图论中的一种算法。由图论可知，图G由两个集合组成：一是顶点的集合V(G)，另一是边的集合E(G)。图G可记做

G=(V，E)聚类开始时，图G为所有的顶点相连的连通图。图中顶点为要分类的化合物，边为多维空间中化合物的距离，当从任一顶点遍历图时，必定将边的集合E(G)分成两个集合T(G)和B(G)，其中，T(G)是遍历图时所通过的边集，B(G)是剩余的边集。显然，G‘=(V，T)是G的子图，称之为连通图的生成树(spanningtree)。所谓最小生成树(minimalspanningtree)即为连通n个顶点的n-1条边的加和为最小的生成树。最小生成树的生成步骤如下：1.设T的初态为空集；2.在连通图上任选一顶点参加到V(T)集合中去。3.将以下步骤重复n-1次：(1)在i属于V(T)，j不属于V(T)的边中选取权值最小的边(i,j);(2)将顶点j参加到V(T)中去；(3)输出i,j及相应权值wij。[例]取本章表11.9中两类的各前10个试样,即试样1-10为多刚毛鸢尾花;11-20为杂色鸢尾花。两个试样间的距离即图G中边的权重以欧氏距离量度。式中，xik为第i个试样的第k个特征，即鸢尾花的花瓣特征。对一组数据，在最小生成树中结点(试样)对及其边的权重值示于表11.19。表11.19最小生成树中结点对及其边的权重图11.11所示为该最小生成树的二维影射图，它由主成分分析(K-L转换)方法所得。假设由结点4和18间(点间距离最大者)断开,那么形成两个局部:左侧全部为类1,即多刚毛鸢尾花;右侧全部为类2,即杂色鸢尾花.由此可见，此最小生成树反映了原数据空间中数据间相邻关系。图11.11不同鸢尾花的最小生成树§11.7显示方法[15-18]如前所述，一种化合物总是定义为多维空间中的一个点，但在多维空间中点的分布(即数据结构)是什么样那么超出了人类视觉判别能力。假假设有方法把这些点投影到低维，如二维或三维空间中那么可借助于图示进行化合物的分类。§11.7.1线性映射(Karhunen-Loeve)主成分分析(PCA)有多种形式，它们构成了多维数据分析的数学根底。Karhunen-Loeve(K-L)转换是一种形式。该种转换是使二维空间中的两个坐标y1和y2(或三维空间中的三个坐标y1,y2和y3)是高维空间中原坐标x1,x2…,xn的线性组合，并且y1和y2(或y1,y2和y3)正交。K-L转换为二维(或三维)坐标中引进均方误差最小的一种转换。如在m维空间中n个点的矩阵：它的协方差阵为：其中为m维空间中n个点的重心将协方差阵对角化，那么得本征矢量矩阵T和本征值λ1,λ2,…,λm其中，λ1≥λ2≥…≥λm;T’=(t1,t2,…,tm)。在如上计算中，亦可用相关矩阵(见后边数字例子)，最终所得结论是一致的。为了在低维，如二维空间进行试样点的显示，取与两个最大的本征值λ1和λ2相对应的本征矢量t(1)和t(2)来计算新的坐标y1和y2：(11.16)但根据情况，有时亦取第1，3；2，3或3，4等与本征值相对应的本征矢量进行计算。运用此种方法所得图形的可信度可用下式判定：其比率越大，说明1和2的代表性越强，即所包含的信息量越大，那么所得映射图接近理想情况。经验说明，其百分比率一般应大于80%。[例]老年性白内障是老年人弱视力的主要致因。在白内障形成时，不仅Na/K的比例有改变，而且其他金属元素也相应有所变化。在作者的实验室[15]测定了白内障晶状体和对照组中Ca,Mg,Na,K,Mn,Fe,Cu和Zn等8种金属元素，其结果示于表11.20中。对原始数据进行了标准化处理。用下式计算相关矩阵：表11.208种金属元素在老年性白内障和正常晶状体中的含量式中rij表示变量xi和xj的相关系数，当i=j时，rij=1。在此例中下脚p=8。运用Jacobi方法计算特征值及与之相应的本征矢量。其本征值为：

λ1=2.3758λ5=0.8863

λ2=1.4998λ6=0.4654

λ3=1.3603λ7=0.1801

λ4=1.1138λ8=0.1185与之相对应的前四个本征矢量为：

t1’=(-0.2329,-0.09168,-0.4712,0.5425,-0.1607,-0.4064,-0.4780,-0.04004)

t2’=(0.5064,-0.2151,-0.3462,0.2501,-0.3585,0.02583,0.4857,0.3871)

t3’=(-0.3975,-0.2944,-0.06075,0.2659,0.3695,0.5641,-0.03175,0.4746)

t4’=(0.04553,-0.8005,0.1987,0.02561,-0.2260,0.1728,-0.06608,-0.8812)利用式(11.16)可计算新坐标y1,y2,…,y8(即主成分PC1,PC2,…PC8)。以主成分1和主成分2张成的二维平面，其投影结果示于图11.12。由此可见，三类(即白内障晶状体，白内障晶状体核和正常晶状体)被清楚地划分在不同的区域。主成分1和主成分2所解释的信息量为47%。假设以主成分1，2和3为坐标，所得三维投影示于图11.13。和图11.12相类同，几种不同的晶体在三维空间中清楚地散布在不同的区间。主成分1，2和3所解释的信息量为65%。图11.12白内障二主成分投影图图11.13白内障三主成分投影图§11.7.2非线性影射(non-linearmapping,NLM)同线性投影一样，非线性投影的目的也是将m维空间的点投影到低维空间如二维或三维空间以使我们易于观察试样的数据结构。即NLM：x→y此处，x和y分别为m维空间和低维，如二维空间中的矢量：

x={xi},i=1,2,…,n(n为试样数)

y={yi},i=1,2,…,n(n为试样数)而设在m维空间矢量xi与xj的距离为及在二维空间中yi与yj的距离为dij=dis(yi,yj)其中，dij’=dis(xi,xj),dij=dis(yi,yj)为xi和xj间的距离表示。在NLM方法中，多维空间中的点经过投影(即在低维如二维或三维空间中)力图保持点与点间的距离不变，在理想情况下dij=dij*，但事实上这是不可能的，因为经过投影之后必然会产生误差：误差函数为采用不同的误差函数E及其距离的测试表达式，那么有各种不同的NLM方法。不管什么样的NLM方法，其操作步骤大致为：步1：设定初始矢量y。从理论上讲，y可以任意给定，不妨将线性投影的结果作为初始矢量y。步2：调解y的结构使得由新y所产生的误差比原来y所产生的误差要小。由现在的数据结构转换到下一数据结构称为迭代。步3：重复步2，当满足如下条件之一时停止迭代：(1)E已到达预先给定值；(2)迭代已到达预先给定次数；(3)操作者对现有数据结构已满意。本节我们主要介绍一下Sammon的NLM算法。令dij*为欧氏距离，即同样，dij亦为欧氏距离，那么其误差函数定义为：在计算中需不断调整在d空间中的n个矢量，直到E到达最小值或预先给定值。也就是说，当E到达预定值时，y即为x在d维空间中的转换矢量。误差函数E是2n个独立变量yij(i=1,2,…,n;j=1,2)的函数，所以在Sammon的算法中，这2n个变量必需同时调整，以得到新的数据结构。为了实现这一点，可以运用最陡下降法或单纯形法来求得E的最小值。假设现在的数据结构为：运用最陡下降法连续计算新矢量这种计算由如下递推公式控制：此处α为校正因子，E为相应于目前数据构形的误差。为yp的一阶偏导数。该组yp’(p=1,2,…,n)将为d空间中新的数据构形。在每一步迭代中，均需计算误差函数对yp的一阶和二阶偏导。其中，一阶偏导为：式中即，C为常数，由于由此可得二阶偏导为：

p=1,2,…,n§11.8综合性数据例子本章我们介绍了模式识别中的各类方法，且对大局部方法已经给出了数字例子。本节我们给出两个比较有代表性的例子，旨在说明在解决实际问题中常常需要综合运用多种方法。另外，为了阅读上的方便，重新给出了计算中的某些公式。[例1]中国茶的模式识别分类[18]从中国不同地区取绿茶，红茶和乌龙茶三类共31个试样，一些化学成分的分析结果示于表11.21。关于茶叶质量的评估由有关专家品尝而得。1．数据的标准化处理将数据阵X(N×M)(N为试样数；M为变量数)的数据元进行标准化处理：其中为均值，sm为标准偏差。表11.21三类茶叶试样化学分析结果C1—C7和H1—H5为绿茶；K1—K4，F1—F7为红茶；T1—T4，S1—S4为乌龙茶，其中C，H，K，F，T和S分别表示中国不同地方之名称的字头。2．系统聚类分析试样间距离采用欧氏距离：类间距离采用Ward算法：上式表示类i与新类(由类j与类k聚合而成)间的距离，其ni表示类i中的试样数。系统聚类结果示于图11.14。图11.14茶叶试样的系统聚类图3．主成分分析对于某一试样i的第k个主成分PC(yik)为其中tmk为相关矩阵(M*N)第k个本征矢量的第m项(相应于第m个变量)。主成分的标准化得分(score)由下式获得：而主成分的装载(loading)的数学表达式为：由装载可以表征原来变量与主成分间的关系。得分和装载的概念可见第五章(即主成分分析和因子分析)。表11.22给出相应于头3个主成分的装载，以及本征值

1---

3等。表11.22茶叶数据主成分分析结果表11.22中公因子方差表示相应于这几个主成分的某一变量所能解释的偏差量，如纤维素可以解释98%，多酚可以解释94%等。对于某一变量公因子方差等于各主成分装载的平方加和。由表11.22可见，头三个本征值能解释的偏差量为96.6%。另外，还可看出，第一个主成分与纤维素类(头三个变量)均为正相关，而与其它三个变量均为负相关。第二个主成分对变量氨基酸最重要，而第三个主成分对变量纤维素最重要。图11.15茶叶三个主成分的投影图以主成分1—3张成三维空间，试样点的投影示于图11.15，而以主成分1和2为坐标的投影示于图11.16。二图均可将绿茶、红茶和乌龙茶三类茶叶清楚地散布在不同的区域。另外，由图11.16似乎可以看出，沿着第一个坐标轴的投影表征了茶叶的质量，即每一类均是沿着坐标1由质量优良的茶叶到质量较差的茶叶。图11.16茶叶二个主成分的投影图4．SIMCA方法分类在SIMCA方法中，对于每一类q构造一主成分模型：假设在类q中有rq个主成分，那么其中Zim(q)意义同前，即分别为均值和标准偏差。对于类q剩余方差为：对于任一试样p，主成分得分ypk由下式计算：其中，Zpm(q)经过标准化处理。此试样矢量由类模型的预测为：而剩余方差为在本例中，按照产地，即C，H，K，F，T和S将三类茶叶分为6个组，对于每组构造一主成分模型。交叉验证的结果说明，每一组的主成分数均为1。图11.17所示为与绿茶类模型的距离图。而距离由前述某试样拟合于类模型时剩余方差(Sp(q)2)公式计算。x轴表示与绿茶C组模型的距离，y轴表示与绿茶H组模型的距离。图中平行于x轴的直线相应于H图11.17距离图“△〞，C组；“〇〞，H组；“●〞，其他组。组模型的临界距离，而平行于y轴的直线相应于C组模型的临界距离。临界距离由F检验测得。判断试样p该否归入某类时的F检验为：其自由度为M-rq和(nq-rq-1)(M-rq)。当计算值F小于某一显著水平下查表得到的F值时，试样应归于某类q。临界距离即为当α=0.05作F检验满足上述条件时的Sp(q)2值。由图可见，所有的试样均得到了正确分类。与绿茶类同，用同样的方法可以对红茶和乌龙茶进行分析。5．质量预测为了定量预测茶叶的质量，我们以主成分(PC1)为横坐标，以质量指数为纵坐标，那么绿茶(C组)和红茶(F组)的回归直线分别示于图11.18和图11.19。图11.18绿茶质量指数回归分析图11.19红茶质量指数回归分析其相关系数分别为0.997和0.983。对于绿茶的H组，红茶的K组和乌龙茶的T组所得结果类似。但对于S组，未得到有意义的结果，因为S2，S3和S4对于主成分1得分几乎相同。[例2]含氟复合化合物中Eu2+跃迁发射的模式识别研究[17]含氟复合化合物中掺入稀土元素是一代新型的发光材料。在作者的实验室曾经进行过含氟复合物中Eu2+发射光谱模式识别的研究。这类化合物的通式为ABmFn，其中A，B分别为某金属原子，F为氟原子。为了模式识别研究，我们提取了如表11.23所示的12个特征。训练集中63个化合物的这些特征的参数及其化合物所属的类列于表11.24。表11.23关于化合物ABmFn的变量表11.24训练集特征续表11.246．变量的选择在模式识别中，变量的选择致关重要，本例中侧重说明变量选择方法的综合运用。(1)方差法变量的方差定义为：此处为试样j的第i变量。很显然，V