GIS空间分析4

第四章 空间分布

空间聚类分析

第一节 概述

聚类分析是一种数字分类方法。对分类对象间的关系用某些相似性度量进行刻划，根据相似性度量进行分类。也有称群分析、簇群分析、簇分析、点群分析、丝分析等，在1975年(苏州)概率统计会上，决定使用聚类分析这一名称。

从数学角度来看，聚类分析是在一个大的对称矩阵中探索诸元素间相关关系的一种数学分类方法。

一．两种不同的模式识别(有无训练集)

有训练集:如判别分析，在已知样本分类的情况下设计判别函数。这些已知分类的样本

称为训练集，也称为有人管理或有教师的分类法。

无训练集:如聚类分析，在设计分类器时，所选用的样本预先不知所属的类别，需要根据样本间的距离或其它相似性的程度来自动地进行分类，也称为无人管理的分类。 二．聚类分析的两种类型

根据分类对象不同，分为Q型(Q-mode)和R型(R-mode)。 Q型

?对样品进行分类。即把不同的物体(如岩石标本、样品等物种或人种)进行比较，?对变量进行分类。即属于同一物种的各种属性，即各类变量(如岩石厚度、岩石

目的是要确定不同物体之间的关系，从而将物体进行归类分群。

R型

成分及各种化验观测数据)进行比较，目的是要不同变量之间的关系，从而对变量进行分类。

三．两种不同的聚类方法(步骤)

根据聚类的方法，分为系统聚类法和分解法两种。

系统聚类法(Hierachical Clustering Methods)是从少到多的聚类法，开始各个样本都各自为一类，以后逐步归并，直至全部样本变成一类。

分解法则相反，开始时全部样本为一类，以后分解，直至各个样本自成一类。 系统聚类法又包括几种不同的方法。这里主要介绍系统聚类法。 聚类分析的结果一般通过聚类图(谱系图、枝状图)反映。

例，在煤田地质勘探过程中，有时煤系含有多个煤层，如果标志层不明显，只用宏观的标志进行煤层对比较为困难，这时就可用聚类分析进行煤层的数字分类，从而达到对比煤层的目的。

做法是在一个煤田或勘探区内，选择若干个煤钻孔，对所有煤层进行工业分析、光谱分析等取得一批实验观测数据(如下表)。

1

指标 化验号 1 2 3 4 5 6 7 8 9 Ag 43.46 37.03 13.90 39.84 36.78 29.73 36.38 38.61 21.72 某煤田煤层化验(部分)成果 煤灰成分 SiO2 52.62 6.10 30.36 71.90 53.58 55.80 67.90 51.06 53.64 Fe2O3 5.70 5.20 33.59 4.00 1.80 7.80 5.40 3.80 10.30 Al2O3 15.24 23.35 28.21 14.59 37.84 23.89 15.84 37.40 25.94 CaO 9.95 2.74 3.09 2.33 3.23 4.19 4.39 3.33 2.72 MgO 5.82 2.86 0.89 1.04 1.09 3.11 2.96 1.42 2.66 SO3 9.60 1.36 1.58 1.55 1.80 2.23 2.19 2.06 1.31 SQg 6.39 7.83 10.13 6.56 5.12 8.51 7.78 6.17 9.95 煤的灰分(Ag)表示出来为：

0 5 10 15 20 25 30 35 40 45 50 Ag(%)

显然：1 点 的Ag>35%较接近，可为一组(群)； 2 点 的Ag在20%与30%之间，可为一组(群)； 3 点 的Ag<15%，自成一组(群)。

如果每个煤层化验两个指标(变量)，如煤的灰分(Ag)和硫分(SQg)，则Ag和SQg为各点的二维坐标，如下图：

SQg (%) 10 9 8 7 6 5

2

4 3 2 1

0 5 10 15 20 25 30 35 40 45 50 Ag (%)

多维空间的点就不能用图形直观的表现出来。表示多维空间点与点之间的疏密关系的量，称为相似性度量(相似性统计量)，如相似系数、相关系数、距离系数、离差平方和增量等。

第二节 数据的规格化处理

为了使观测数据以等同的量级出现，必须规格化处理，通常的方法有：标准化、正规化、均值化及对数变换等。

１．数据标准化

设有n个样品，每个样品测量了m项指标(变量)，得出如下原始数据矩阵，为

X??xij??x11?x21??????xn1x12x22?xn2?x1m??x2m???????xnm?

i?1,2,?,n（样品个数）j?1,2,?,m（变量个数）xij表示第i个样品、第j个变量的观测值。

设变换后的数据记为zij，则 zij?xij?xjsj(i?1,2,?,n,j?1,2,?,m)

1n 其中xj??xij,sj?ni?1 写成矩阵形式为

1n(xij?xj)2 ?n?1i?1 Z?zij???z11z12?zz2221???????zn1zn2?z1m??z2m???????znm?

则称Z?zij为标准化数据。

?? 3

注：若所取样品构成的变量服从正态分布，则标准化后的数据Zij~N(0,1)。 ２．数据正规化 所谓 数据正规化，就是通过极差变换，把原始数据矩阵中的任何一列的最小值化为0，最大值化为1，其余介于0与1之间。记 Zij??xj??min?xj?max?xj?xij?min(1?j?m)

minxj分别为第j个变量的最大值与最小值。 其中maxxj、写成矩阵形式为：

???? Z?zij???z11z12?zz2221???????zn1zn2?z1m??z2m??

?????znm?

第三节 相似性度量

能够度量变量(或样品)之间相似性程度的数量指标，称为相似性度量(统计量)。常用的有：相关系数、相似系数、距离函数、误差(距离)平方和增量等。

一、 相关系数

在R型分析中，定义第j个变量和第l个变量之间的相关系数为：

?jl?Cov(xj,xl)?j?l (理论相关系数)

Cov为协方差，?j、?l为标准差。 实际应用中，定义相关系数为：

rjl?i?1n?(xnij?xj)(xil?xl)n

i?1?(xijij2?xj)?(xil?xl)2i?1其中，xj1?ni?1?xn,xl1?ni?1?xnil(j,l?1,2,?,m)

相关系数矩阵为：

4

?r11?rR??21????rm1?r1m?r22?r2m?? ?????rm2?rmm?r12为一m×m阶矩阵。

在Q型分析中，定义第i个样品和第k个样品之间的相关系数为：

rik?j?1m?(xmij?xi)(xkj?xk)m

j?1?(xmj?1ij2?xi)?(xkj?xk)2j?1其中，xi1?mj?1?xr12r22?rn2mij,xk1?m?xkj(i,k?1,2,?,n)

相关系数矩阵为：

?r11?r21R??????rn1?r1n??r2n??

?????rnn?为一n×n阶矩阵。

?1?r?1。r大?关系密切。

二、相似系数(cos?)

Q型：

Q型中相似系数(cos?)用于样品间的比较，夹角?愈小，表示两个样品愈相似；夹角?愈大，表示两个样品愈疏远。下面讨论m=2情形。

若有n个样品，每个样品测得两个变量(x1,x2)，则每个样品可用在x1ox2坐标系中的点表示出来。若任取A、B两个样品，则A、B在x1ox2坐标系中用相应的A、B两点表示出来，如下图所示：

x2

A(x1A,x2A) B(x1B,x2B) a b

?1 ?2

5

O x1

二维样品空间图

OB的模长分别为a、b，夹角为?，则： 设OA、 co?s?cos?(s1co?s2?sin?1sin?2 1??2)?co? sin?1?x2Axxx,sin?2?2B,co?s1?1A,co?s1?1B, abaaxxxx∴ cos??1A1B?2A2B?a?ba?bj?1?x2jAxjBj?1?x?x2jAj?1222jB

推广到m维空间，即每个样品测得m个变量(x1,x2,?,xm)，xi(xi1,xi2,?,xim)、

xk(xk1,xk2,?,xkm)为m维空间中的两个向量，由于(xi,xk)?|xi||xk|cox?ik (内积)

∴ cox?ik?(xi,xk)|xi||xk|

据内积定义(xi,xk)?代入(3)式中，得 cos?ik?令 sik?cos?ik?xj?1mijxkj, |xi|=

?xj?1m2ij,| xk|=

?xj?1m2kj

?xj?1mijxkj?x?x2ijj?1j?1mm2kj i,k?1,2,...,n

(i,k?1,2,?,n)

?由于cos?ik为向量xi和xk的夹角，在标准化条件下，其取值范围是0~360?，因此相似

??系数(cos?ik)的取值范围为－1—1之间；在正规化条件下，其取值范围是0~90，因此相似系数(cos?ik)的取值范围为0—1之间。

显然，sik大，表示i,k两个样品相似系数愈大。 对n个样品可得

?s11s12?ss2221?S?(s)? ik?????sn1sn2为一n×n阶矩阵。

?s1n??s2n??

?????snn? 6

R型：

R型中相似系数(cos?)用于变量间的比较。

cos?jl??xijxili?1n?x?x2iji?1i?1nn2il j,l?1,2,...,m

此时组成一个m×m阶矩阵。

,sj?sl?1，则有： 在标准化数据下，由于xj?xl?0 rjl?i?1?x2ilnijxilni?1?xn?cos?jl2il

i?1?x即R型标准化数据的相关系数等于相似系数。 三、 欧氏距离

用于Q型。

欧氏距离：在直角坐标系中，两个样品相互之间的距离。

若两个样品A、B各自只有1个观测值xA、xB时，dAB?xA?xB 每个样品有x1、x2两个变量，则 x2

dAB???22(?x1)?(?x2)22(x1A?x1B)?(x2A?x2) A

j?1?(x2jA2?xjB) 推广到m维空间中，即设有n个 B 样品，每个样品有m个变量(x1、x2 ?x1 ?x2 、? 、x1)，若有i、k两个样品点， x1 其坐标分别为：

?xi(xi1,xi2,?,xim) ??xk(xk1,xk2,?,xkm)则i，k两点间的距离为：

dik??(xj?1mij?xkj)2 i,k?1,2,...,n

7

D?d11?d21?（dik)??????dn1d12d22?dn2?d1n??d2n??

?????dnn?

四、斜交距离(系数)矩阵

用于Q型。

斜交距离实际上是斜交(仿射)坐标系下两点之间在欧氏空间的距离，或说是校正到欧氏空间的距离，是由斜交分量表示的欧氏距离。

odik?j?1??(xl?1mmij?xkj)(xil?xkl)cos?jl

nncos?jl?

i?1?xnijxili?1?x?x2iji?12il

五、离差平方和增量(?E)

离差平方和增量是以组内方差为基础推导的，多用于Q型聚类分析。 定义 xijk (k=1,2,…,G;j=1,2,…,m;i=1,2,…,nk) 为第k类中第i个样品的第j个变量的计量；

xjk1?nki?1?xnknkijk

为第k 类中关于变量j的平均值 1. Ek?j?1i?1??(xmijk2?xjk)

称为第k类的组内离差平方和。

例： 有5个样品，每个样品两个变量，分成三类，数据如下：

??3??1??2?4??5?? -8 A??? 5 B? 4?? 1 C?-2??nA=1, nB=2, nC=2, EA=0, EB=1, EC=5。

2. 在系统聚类法中，开始时每个样品自成一类(每类包含一个样品)，此时对每类都有

xjk?xijk

8

Ek=0

现在考虑合并问题，假如原有两类p、q要合并成新类t,各类样品数分别为np, nq, nt,平均值分别为xjp,xjq,xjt,(j=1,2,…,m)，显然nt=np+nq,称

?Epq?Et?(Ep?Eq)

为离差平方和增量。

如上例中B与C合并，先算出EBC=40, 则?EBC?40－1－5＝34 由

Et?mj?1i?1np??(xmntijt?xjt)?m2j?1i?1np??xmnt2ijt?nt?x2jt

j?1mmEp?j?1i?1mnq??(xijp?xjp)?2j?1i?1mnq??x??x?2ijp?np?x2jp

j?1mEq?j?1i?1??(xijq?xjq)?mnt22ijqj?1i?1?nq?x2jq

j?12ijp和

j?1i?1??xm2jp2ijtj?1?(?xi?1mnp?i?1?xnq2ijq)

可得： ?Epq??ntj?1?xm2jt?np?xj?12?nq?xjq

j?1m3. 上式进一步整理化简得：

?Epq?npnqnp?nq?(xj?1mjp?xjq)2

该式说明p、q两类合并后离差平方和增量?Epq与两类中的各变量的重心(平均值)的距离之差成正比，并以

npnqnp?nq为比例系数。

4. 当p、q两类合并成t时，其它类(设为r)与p、q的离差平方和增量?Erp、?Erq要刷新为

?Etr,可用刷新方程 ?Etr?1[(nr?np)?Erp?(nr?nq)?Erq?nr?Epq]

nt?nr5. 离差平方和增量法聚类过程

设样品数为n，每个样品有m个变量，开始时每个样品自成一类，计算初始离差平方和

9

增量矩阵(?Eik?Eik?Ei?Ek?Eik1?2j?1?(xij?xkj)2?m12dik)： 2?E(0)??E11??E21???????En1?E12?E22??En2??E1n???E2n??

??????Enn?从中选择?Epq最小的两类归并为一类(t)，再求类t与其它类的误差平方和增量(利用刷新方程)，?E(0)(n×n阶)??E(1)((n-1)×(n-1)阶)。再从?E(1)中选择?Epq最小的两类归并，直至全部样品归并成一类为止。要注意的是每归并一次?E的阶就减少1，到最后变成2×2阶。

例： 有一批样品，原始数据为

8??3－?15????24? ??41???2??5－?

初始离差平方和矩阵

1 2 3 4 5?173145?04120??22?2565?01??22? ? ?1345??0?22??05???0??(0) ?E(0)可以看出，其中?E23=1为最小，表明类2与3最相似，将它们合并成新类，

1 2_3 4 5?317?04120??3?37109?

0?? ?33??05???0??? ?E(1)现在是?E45＝20最小，故将类4与5合并成新类，

10

(1)

1 2_3 4_5?317?039?? 3? ?034???0???? ?E(2)再将类2_3与4_5合并成新类，

1 2_3_4_5?080? ?E(3)? ?0???

至此，只剩下2类，1于2_3_4_5，它们的离差平方和增量为80,故不需计算，将它们合并成

一类即可。

上述过程可以写成联接表。

?E L(联接序号) 联接样品 1 1 2～3 2 5 4～5 3 34 2_3～4_5 4 80 1～2_3_4_5 最后画出聚类图

第四节 系统聚类法(Hierachical Clusting Metjods)

一、一般步骤

系统聚类法首先认为n个样本(样品或变量)自成一类，此时各样本之间的相似性度量为初始相似性度量矩阵；然后选择最相似的两类合并成一个新类，计算新类和其它类的相似性度量，再将最相似的两类并成一类，这样每次减少一类(相似性度量矩阵缩小一阶)，直至所有的样品合并成一类为止。 它的一般步骤为：

1．开始时每个样本自成一类，共有n类，用数1,2,…,n分别标记这n类。

2. 从相似性度量矩阵中找出最相似的两类，分别用p和q(p

3. 将p、q合并，新类记为p(保留p行p列)，删除矩阵中属于q的行与列(q行q列)，标记p的合并结果并刷新(更新)相似性度量矩阵中类p和所有其余现存类之间的相似性度量。

4. 重复第2步和第3步，总共进行n-1次，此时全部成员都包括在一个类内。每次将合并的类和联接它们的相似性度量记录下来，得到聚类分析的一份完整记录。

各种不同的系统聚类法的区别无非在于第3步中刷新相似性度量矩阵的方法而已。所谓

11

最相似的两类，对于距离类度量，是距离最小的两类；对于相关系数和相似系数，是相关系数或相似系数最大的两类。 二、系统聚类法各论

1. 最短距离法(近邻连接法)

以两类中最相似的两个样本的相似性度量作为该两类的相似性度量。 一般公式：

1) 对欧氏(斜交)距离

str?min(spr,sqr)

2) 对相关系数和相似系数

str?max(spr,sqr)

递推公式：

1) 对欧氏(斜交)距离

str? 2) 对相关系数和相似系数

111spr?sqr?spr?sqr 222111spr?sqr?spr?sqr 222str?2. 最长距离法(远邻连接法)

以两类中最不相似的两个样本的相似性度量作为该两类的相似性度量。 一般公式：

1) 对欧氏(斜交)距离

str?max(spr,sqr)

2) 对相关系数和相似系数

str?min(spr,sqr)

递推公式：

1) 对欧氏(斜交)距离

str? 2) 对相关系数和相似系数

111spr?sqr＋spr?sqr 222111spr?sqr－spr?sqr 222str?3. 类平均法

两类之间的相似性度量取为该两类中不同类任意两个样本的相似性度量的平均。 如类p和类q的相似性度量

spq?1npnq1?i?np1?k?nq?sx(p),x(q)

ik其中xi(p)为类p中的样本，xk(q)为类q中的样本。

12

x(q) 1x(p) 1x(p) 2x(p) 3 递推公式：

x(q) 2x(q) 3x(q) 4

npnqspr?sqr

np?nqnp?nqstr?4. 重心法

在聚类分析的每一阶段将平均数(或重心)最接近(相似)的两类合并，称为重心法。

spq?sxp,xq 等于两类均值之间的相似性度量。

若p、q合并为新类t，则新类t的重心

wt?npwp?nqwqnp?nq 或 xt?npxp?nqxqnp?nq

类t与其它类r的相似性度量 str?sxt,xr

如果相似性度量是相关系数和相似系数，目前重心法还未有递推公式；所以每次合并后都要重新计算新类t与其它类的相似性度量。 对于欧氏(斜交)距离的平方，有递推公式：

2令str?dtr，spr?dpr，sqr?dqr,spq?dpq，则

222str?npnp?nqspr?nqnp?nqsqr-npnq(np?nq)2spq

重心法的缺陷是可能发生“反转”现象，即在聚类过程中，随着类的扩大，相似性反而增强，这是不合理的。其它几种方法无此现象，即在聚类过程中，相似性都是单调下降的。 5. 离差平方和增量法 见上。 三、系统聚类法小结

这五种系统聚类法都是常用的，它们的共同点是每归并一次，都要改变相似性度量矩阵，除重心法相关系数和相似系数外，都可用递推公式。这些递推公式可写成一个统一的格式： 设p和q归并成类t，则类t与其余任一类r的递推公式为：

str??pspr??qsqr??spq??spr-sqr

式中参数见下表。

系统聚类法递推公式参数表 方法 ?p ?q ? ? 距离类 相似系数类 13

最短距离法 (近邻连接法) 最长距离法 (远邻连接法) 类平均法 重心法* 离差平方和 增量法 1 21 2npntnp 1 21 2nqntnq 0 0 0 －1 21 2 0 1 21－ 2 0 / / ntnp?nrnt?nr ntnq?nrnt?nr nnr?nt?nr ?npnq2t 0 / *重心法仅适用于欧氏(斜交)距离的平方，其中str、spr、sqr、spq均为距离平方。 一般说来，最短距离法适用与条形的甚至S性的类，其它方法适用于椭球形的类，在使用时应根据实际问题加以选择。

14

第五节 聚类分析应用实例

例1 某地煤层中采了10块煤样，作了显微组分百分含量分析，6个变量分别代表镜质体、壳质体、半丝质体和丝质体、其它惰性体、黄铁矿、其它矿物。 97.0 1.7 0.4 0.2 0.5 0.2 96.5 1.5 0.3 0.2 1.0 0.6 78.8 10.4 5.0 2.4 1.3 2.1 80.9 7.9 4.3 2.2 2.9 1.8 66.3 15.1 9.1 3.4 2.1 4.0 65.0 12.1 8.2 3.1 5.6 6.0 35.9 14.5 17.4 11.6 0.7 19.9 50.6 15.8 12.2 4.0 4.6 12.8 10.2 2.0 84.5 1.2 1.4 0.3 16.0 0.5 82.7 0.0 0.3 0.5 1. Q型分析

1) 原始数据正规化，采用相似系数，最长距离法 原始相似系数矩阵 1

0.9953 1

0.7967 0.8018 1

0.8007 0.8332 0.9488 1

0.5862 0.6008 0.955 0.901 1

0.4992 0.5574 0.8133 0.9117 0.8794 1

0.2243 0.2318 0.6291 0.534 0.7377 0.6094 1

0.3715 0.418 0.7896 0.818 0.9052 0.9472 0.7927 1

0.0179 0.0344 0.1642 0.1933 0.2315 0.2683 0.2346 0.2783 1

0.069 0.0677 0.1042 0.0957 0.1204 0.0965 0.1382 0.1173 0.9674 1 Q型相似系数最大距离法联接表 L(联接序号) 联接样品 相似系数 1 1～2 0.995 2 9～10 0.967 3 3～5 0.955 4 6～8 0.947 5 3～4 0.901 6 3～6 0.790 7 3～7 0.534 8 1～3 0.224 9 1～9 0.018

15

Q型相似系数最长距离法聚类图

2) 原始数据标准化，采用斜交距离，重心法 原始斜交距离矩阵 0

0.2781 0 2.472 2.31 0

2.577 2.352 0.7798 0

4.009 3.835 1.543 1.681 0 5.08 4.846 2.88 2.51 1.762 0

8.489 8.388 6.241 6.582 4.998 5.423 0

6.381 6.184 3.963 3.893 2.457 1.825 3.787 0 5.15 5.159 5.16 5.334 5.664 6.522 8.322 7.092 0 5.053 5.105 5.466 5.677 6.17 7.139 8.965 7.791 0.9303 Q型斜交距离重心法联接表 L(联接序号) 联接样品 斜交距离 1 1～2 0.278 2 3～4 0.780 3 9～10 0.930 4 3～5 1.566 5 6～8 1.825 6 3～6 2.861 7 1～3 3.962 8 1～9 5.497 9 1～7 6.328

Q型斜交距离重心法聚类图

16

0

3) 原始数据标准化，采用离差平方和增量法 原始离差平方和增量矩阵 0

0.041 0

1.506 1.442 0

1.727 1.407 0.4801 0

3.841 3.654 0.5571 1.005 0

6.715 5.954 3.188 1.729 2.055 0

14.33 14.19 8.648 10.07 6.669 9.629 0 8.80 8.142 3.931 3.299 2.016 1.017 5.58 0

7.408 7.273 6.421 6.481 6.954 8.855 13.69 9.478 0

6.618 6.647 6.577 7.049 7.739 10.68 14.95 11.26 0.2973 Q型离差平方和增量法联接表 L(联接序号) 联接样品 离差平方和增量 1 1～2 0.041 2 9～10 0.297 3 3～4 0.480 4 3～5 0.882 5 6～8 1.017 6 1～3 4.861 7 6～7 9.801 8 1～6 16.62 9 1～9 20.59

2. R型分析

1) 原始数据不变换，采用相关系数，类平均法 原始相关系数矩阵 1

0.0214 1

-0.8794 -0.4621 1

-0.239 0.7083 -0.2046 1

0.0875 0.5714 -0.3045 0.0866 1

-0.2522 0.7303 -0.2108 0.9338 0.2415 1

17

0

R型相关系数类平均法联接表 L(联接序号) 联接样品 1 4～6 2 2～4 3 2～5 4 1～2 5 1～3

相关系数 0.934 0.719 0.300 -0.096 -0.412

R型相关系数类平均法聚类图

例： 现有一个对样品的聚类分析, 以欧氏距离作为相似性度量，采用最短距离法进行聚类，已进行了二次归并, 联接表如下：(1,??3,4,?为样品号)

L ????? D pq 聚类样品 1 1?4 ??? 2 ??5 ? 3 4 1,? ? ???

?054.5??07? D =

???0???

请完成该聚类分析,并画出聚类图(要求列出算式)。 解：第3次归并为1,? 与3,?，由最短距离法递推公式

111Str?Spr?Sqr?Spr?Sqr

222得：

?,???,?,???2 D =??05?? 0??

故联接表为： L 聚类样品 1 1～4 2 ?～5 3 ??4～3,? 4 1,?,???～2 聚类图：

??D pq ? ? 4?? ? 18

1 ? ? ? 2

19

趋势面分析

第一节 趋势面分析的概念及类型

一． 趋势面分析的概念

趋势面分析是拟合空间曲面(或超曲面)的一种统计方法，是回归分析的直接推广。在地质上，趋势面分析可将地质变量(特征)区分为区域性变化分量和局部性变化分量，从而研究地质变量(特征)的空间分布和变化规律。 1. 趋势：指事物发展或变化的主体，受总的规律的支配。在地质上往往由区域构造、区域岩相等大区域因素所决定。

2. 局部异常：受局部因素支配，为事物发展的特殊部分。在地质上往往是注意的重点。以成矿为例，元素的平均值是很低的，唯有成矿(富集)才能形成可开采的矿。这些富集点往往被看作局部异常点。

当然，趋势与局部异常是相对而言的，视研究区、研究对象不同相互转化。 3. 随机性变化：由随机因素形成的偏差(白噪声)，包括取样和分析误差。 地质变量的观测值，如煤层厚度、矿石品位等地质变量?x地?，可表示为

x地?x区?x偏x偏?x局?x随

?x区—趋势部分??x偏—偏差部分其中 ? ，即 x地?x区?x局?x随

?x局—局部偏差?x—随机偏差?随例1． 煤层厚度变化趋势图

逐渐变薄

20

例2． 构造引起煤层局部变化图

煤厚在小断层处变厚

4. 趋势面: 反映地质变量在空间趋势变化的曲面。一般指二元函数，如煤层厚度、煤层顶底板标高、构造面高程等，这是一个真正的曲面。但也有指三元函数的，此时称为超曲面。 5. 趋势面分析：对于空间一组观察值，排除随机性干扰，找出反映主要变化规律的趋势

?(面)方程z?f(x,y)，求出局部异常。利用得出的结果结合地质条件，分析趋势变化和局部异常的地质意义，从而指导进一步的地质工作(勘探或找矿)或对地质过程作出合理解释。 趋势面分析一般作两种图，趋势图和偏差图(剩余图)(ei?zi?zi)。直接用原始数据插值得到的称为实测图，实测图减趋势图就得到偏差图。 二． 趋势面分析的类型

1． 据用于拟合趋势面的函数划分： (1) 多项式趋势分析(常用)

(2) 调和趋势分析(Fourier趋势面分析)

(3) 其它函数趋势分析(指数、对数、幂函数等) 2． 多项式趋势面分析的进一步划分：

根据多项式函数中自变量的个数，可分为一元、二元、三元趋势分析，每一种又可分为一次、二次、三次、… 等。

?(1) 一元趋势分析 z?f (x)反映平面上曲线的变化趋势。 一次：z?b0?b1x二次：z?b0?b1x?b2x

???

2

3三次：z?b0?b1x?b2x＋b3x?2

21

………………………………………………

?(2) 二元趋势分析 z?f(x,y)一次：z?b0?b1x?b2y二次：z?b0?b1x?b2y＋b3x三次：

反映三维空间曲面变化趋势，是“真正的”趋势面分析。

?

2??b4xy?b5y2

?z?b0?b1x?b2y＋b3x2?b4xy?b5y2 ＋b6x?b7xy?b8xy?b9y3223………………………………………………

?(3) 三元趋势分析 z?f(x,y,w)一次：z?b0?b1x?b2y＋b3w二次：

??z?b0?b1x?b2y?b3w＋b4x2?b5xy?b6xw?b7yw?b8y2?b9w2三次：

?z?b0?b1x?b2y?b3w＋b4x2?b5xy?b6xw?b7yw?b8y2?b9w2＋b10x3?b11y3?b12w3?b13x2y?b14xy2?b15y2w?b16yw2?b17xw2?b18x2w?b19xyw………………………………………………

一、二、三元趋势方程的次数、自变量个数(不含b0)与极值数关系表 自 变 量 数 次数 极值最大曲线 曲面 超曲面 数 1 1 2 3 0 2 2 5 9 1 3 3 9 19 2 4 4 14 34 3 5 5 20 55 4 6 6 27 83 5 … … … … … 2p p p(p+3)/2 p(p+6p+11)/6 p-1

22

第二节 多项式趋势面分析的数学模型及计算

一． 二元一次趋势面方程及计算

设某一地质变量，在平面上有n个观测值zi与地理横坐标(y)、纵坐标(x)具有相关关系，其表达式为：

? z?f(x,y)

(zi,xi,yi),i?1,2,?,n构成趋势面分析的基本数据，一次趋势面方程为

?z?b0?b1x?b2y各观察点的趋势值为：

zi?b0?b1xi?b2yi据最小二乘原理，要求所配曲面与实测数据的离差平方和

n?22 Q??(zi?zi)??[(zi?b0?b1xi?b2yi)]?mini?1i?1n(1)

?

(2)

n??Q??2?(zi?b0?b1xi?b2yi)?0?i?1??b0n???Q∴???2?(zi?b0?b1xi?b2yi)xi?0?bi?1?1n??Q??2?(zi?b0?b1xi?b2yi)yi?0??i?1??b2(3)

23

(3)式经整理得正规方程组，为：

nnn??b0n?b1?xi?b2?yi??zii?1i?1i?1?nnn?n2 ?b0?xi?b1?xi?b2?xiyi??xizii?1i?1i?1?i?1nnn?n2?b0?yi?b1?xiyi?b2?yi??yizii?1i?1i?1?i?1(4)

矩阵形式为：

??n?n ??xi?i?1?n??yi??i?1?x?xi?1nii?1i?1nni2i?xyi??b0??n?yz???i???i?i?1????ni?1?nb1?xiyi?????zixi??????i?1?i?1n??b??n?22yi?????ziyi????i?1???i?1?????n(5)

系数矩阵A为对称、正定。解(5)得b0、b1、b2，从而得二元一次趋势面方程并可计算出各观察点的趋势值：

zi?b0?b1xi?b2yi?(i?1,2?,n)(6)

它就表示空间一个平面，其等值线图为一组平行线。 偏差值 ei?zi?zi?(i?1,2,?,n)

二． 二元二次趋势面分析

二元二次趋势面方程为：

?z?b0?b1x?b2y?b3x2?b4xy?b5y2各观察点的趋势值为：

(7)

?zi?b0?b1xi?b2yi?b3xi2?b4xiyi?b5yi2

i?1,2,...,n 24

?Q??(zi?zi)2若

i?1nn

22??[(zi?b0?b1xi?b2yi?b3xi?b4xiyi?b5yi)]2?mini?1(8)

?Q?0?bj(j?0,1,2,3,4,5)(9)

整理得六元一次方程组

??n?n?xi??i?1?n??yi?i?1?n2??xi?in?1?xyii??i?1?n??yi2??i?1?n???zi?i?1?xi?xi?1ni?1ni?1n2in?yi?xyii?1ni?1nin?xii?1nn2?xi?1n?xiyi?xi?1ni?1ni?1n2in3iyi2?xiyi?xiyi2i?1n?yii?1n2i2?xiyi2i?1n2?xiyi?xi?1n3i2?x3i?xi?1n2i?1nyi?xyii?1nii?1i?1n?xiyi?yi?1ii3i?xiyi3i?1n?x2i4iyi2in?xi?1niii?1i?1n?xiyii?1nyi2?xyii?1iii3in??b???0?i?1???n2xiyi??b1?????i?1n???3yi??b2?????i?1n22???xy?ii?b3

??i?1?n??3?xiyi?b4?????i?1???n4yi??b5???i?1????yin2?zx?zy?zxii?12?zxyi?1n2?zy?ii?i?1?T(10)?、b5，得趋势面方程并可计算出各观察点的趋势值为： 解得b0、b1、?22zi?b0?b1xi?b2yi?b3xi?b4xiyi?b5yi(i?1,2?,n)(11)

三． 二元p次趋势面分析

对于高次趋势面方程，可借助于类似九九乘法表的形式，得出它们的正规方程组的系数矩阵

25

1 x y x xy y x xy xy y ……

2232 23

1 ?1 ?xi ?yi ?xi ?xiyi ?yi ?xi ?xiyi ?xyi ?yi

2

x ?xi ?xi y

2

x xy

2

y

3

x

2

xy

2

xy

3 3342345332456

y?yi ?xiyi ?yi ?xiyi ?xiyi ?yi ?xiyi ?xiyi ?xyi ?yi . . .

与回归分析一样，可完全用矩阵算法求解，设

2

2

3

2

2

3

?1??1X＝?...??...?1?p?...y12p?x2y2x2xyy...y22222?? .....................?.....................?p?xnynx2xnyny2...ynnn?n?[p?(p?3)/2?1]x1y12x1x1y12y1?b0?z?1??b??z??1?2Z??? B??b2?

?...???...?????zn?n?1??bk??[p?(p?3)/2?1]?1则解为

B?(XX)XZ

实际上，这是矛盾方程组 XB＝Z的最小二乘解。 二元趋势面的计算方法不难推广到三元的情形。

T?1T 26

第三节 趋势面分析的统计度量和限制

一、趋势面系数的假设检验

其实质与回归系数的假设检验完全一致。 １．全体系数检验

趋势面方程的全体系数检验，是检验地质变量(z)与自变量及其高次幂和交叉乘积幂之间线性相关是否显著。

同样有： ST?SR?SD 总离差平方和 ST?i?1n?(zni2 ?z)回归平方和 SR?i?1n?(z?i?z)2

偏差平方和 SD?i?1?(zi?2 ?zi)设p次趋势面分析的系数个数为k(不包括b0)，

零假设(因变量和自变量及其高次幂和交叉乘积幂之间不相关，则偏回归系数全为0)：

H0:b1?b2???bk?0

统计量：

F?SRk~F(k,n?k?1)

SD(n?k?1)给定?，若F?F?，趋势方程面不显著；反之若F?F?，则认为因变量与自变量之间线性相关，趋势方程面显著。 ２．趋势面方程的逐次检验

即趋势面方程增高一次的偏回归系数检验。

设p+1次的系数个数为k, p次的系数个数为q，可分三次检验： (1) 对p+1次k个系数进行检验(全体系数)

27

F(p?1)S(p?1)Rk?(p?1)~F(k,n?k?1) Sn?k?1)D((2) 对p次q个系数进行检验(全体系数)

F(p)S(p)Rq?(p)~F(q,n?q?1) SD(n?q?1)(3) 对p次到p+1次增加系数进行检验

H0:bq?1?bq?2???bk?1?bk?0

统计量：

P?1)P)(S(?S((k?q)RR)F?~F(k?q,n?k?1) P?1)S((n?k?1)D给定?，若F?F?，接收H0，表明趋势面方程增加一次，回归效果不显著增加；若F?F?，则否认H0，表明趋势面方程增加一次，回归效果显著增加。 ３．趋势面方程中单项偏回归系数检验

又分两种情况(设p次趋势面方程有k个系数)。 (1) p次趋势面方程中去掉第j项后的显著性检验 先求出S(p)(p)R，k?1，D，k?1

SH0:b1?b2??bi?...?bk?0 (i?j) S(p)R,k?1(k?1)F?(p)~F(k?1,n?1)

SD,k?1(n?k)(2) p次趋势面方程中第j项偏回归系数的显著性检验 先求出SR,(j)?S(p)(p)R?S(p)R，k?1，S(p)D，k?1

H0:bj?0

S(p)R,(j)1F?(p)~F(1,n?k?1)

SD(n?k?1)

28

对系数的单项检验，类似于逐步回归的情况。所以，我们可以把趋势面分析与逐步回归分析结合起来，或者干脆用逐步回归的方法(此时系数矩阵要作一些改变，即去掉b0)，挑选真正对描述趋势变化有用的项进入方程。

二、拟合度

不同次数的趋势面对原始资料逼近程度显然是不一样的。趋势面对原始资料的逼近情况，在数学上可用离差平方和的变化来表示 Nc?SR?100% ST 当Nc=100%,则趋势值在所有数据点上和观察值完全吻合，这种拟合很少出现。但当数据点数与趋势函数的项数(含b0)相同时，就会得到完全的拟合。这样的拟合是没有什么必要的，因为趋势面分析的一个重要方面就是求偏差值。

拟合度只是一个参考指标，这是因为拟合度只与观察点的趋势值及偏差有关。至于观察点分布是否合理，是否控制了测区，那是另一回事，拟合度不能反映。 拟合度一般取40～60％即可。

第四节 趋势面分析在地质上的应用实例分析

例1． 在15?15km范围内，选择15个观测点测量某沉淀物的粒径，其数据如下表(坐标

单位：千米；粒径单位：毫米)。根据表中数据，求出一次趋势面方程，并计算各点的一次趋势值和偏差值。

点1 2 3 4 5 6 7 8 9 10 号 Xi 2 2 2 4 5 5 7 7 8 10 Yi 3 10 13 1 8 14 3 6 11 8 zi 1.9 2.3 1.1 2.6 2.2 1.8 3.5 3.1 1.3 1.2 解：先求出一次趋势面方程 z?b0?b1x?b2y ?b0?3.10,b1?0.054,b2??0.101 z(x,y) = 3.10078 -0.0538957x -0.100554y

∴

11 12 13 14 15 211 12 12 12 15 13 3 6 10 13 1.4 1.7 1.8 1.2 1.0 ? 29

?zi?3.10?0.054xi?0.101yiST?SR?i?1n?（z?i?1ni2?z)?7.83,fT?15?1?14

?（zi2?z)?3.75,fR?2

Nc?SRST?48.72%SD?ST?SR?7.83?3.75?4.08,fD?15?1?2?12F?SRk3.752??5.51SD(n?k?1)4.0812给定??0.01,F0.01(2,12)?6.93,??0.05,F0.05(2,12)?3.89 ∵ F0.05(2,12)?F?F0.01(2,12)

∴ 趋势面方程在0.05水平下显著，但在0.01水平下不显著。 再求二次趋势面方程 z?b0?b1x?b2y＋b3x＋b4xy＋b5y z(x,y) = 2.50336 + 0.166575x -0.112533y -0.0161639x2+ 0.00332764xy+0.000124488y2

SR=4.48 SD=3.35

?22Nc?SRST?57.22%

SRk4.485F???2.41

SD(n?k?1)3.359

??0.05,F0.05(5,9)?3.48

5,9) ∵ F?F0.05(∴ 趋势面方程在0.05水平下不显著。

14.0012.0010.008.006.004.002.000.000.002.004.006.008.0010.0012.0014.00 实测图(曲面样条函数插值)

30

14.0012.0010.008.006.004.002.000.000.002.004.006.008.0010.0012.0014.00 一次趋势面图

14.0012.0010.008.006.004.002.000.000.002.004.006.008.0010.0012.0014.00 二次趋势面图

14.0012.001.201.000.8010.000.600.400.200.00-0.20-0.40-0.60-0.80-1.00-1.20-1.40-1.60-1.808.006.004.002.000.000.002.004.006.008.0010.0012.0014.00

一次趋势面偏差图(剩余图)

31

14.0012.0010.008.006.004.002.000.000.002.004.006.008.0010.0012.0014.00二次趋势面偏差图(剩余图)

32

-0.40-0.60-0.80-1.00-1.20-1.40-1.60-1.80-2.00-2.20-2.40-2.60-2.80-3.00-3.20-3.40-3.60-3.80-4.00-4.20-4.40

第四章 判别分析(Discrimnant Analysis)

第一节 判别准则

一．判别分析的概念

自然界中分类包括两个方面的内容：其一是研究对象存在着几种类型，即能分为多少类；其二，在研究对象类型数目已知情况下，某一研究个体应该属于哪一类。第一种情况属聚类分析，第二中情况则属判别分析。

判别分析的基本思想：是将研究对象(某一个体)的各种地质特征，同它可能归属的各个类型的地质特征进行对比，以决定其应该归入哪一类。

应用判别分析必须以事先知道存在几个母体为前提，参加建立判别分析的样本必须知道其归属；而聚类分析则不需要这些条件。所以称聚类分析为无训练样本的判别，而判别分析则称为存在训练集条件下的分类。

判别分析的关键在于建立判别函数，这样的判别函数应能有效地区分两类(或多类)事物。拿直观的图来解释，就是使各类重叠区域尽可能地小。

二．判别函数

33

判别函数的形式是一个或几个变量的线性组合(称为线性判别函数)

y?c1x1?c2x2???cpxp?j?1?cxjpj (1)

这样的一个线性组合，比起单个变量来，往往能更好地分辨事物的种类。如上图有A、B两个总体，从中抽取两组样品，每个样品有两个变量，两类同一变量之间，总有些重叠部分。利用两个变量的线性组合构成一个判别函数后，其重叠部分比x1与 x2单个变量的重叠部分都要小。

三．两种判别准则

判别函数的建立有不同的准则: 费歇(Fisher)准则、贝叶斯(Bayes)准则、最小二乘(Lms)准则、库巴克(Kullback)准则、不确定性准则等。本章介绍前二种。

Bayes准则是将m维欧氏空间R(样品是这空间的一个点)划分成G个互不相交的子空间。这样任一个体就可知道它的归属。任何一种划分都可能存在错分，使错分概率最小的分法就叫做Bayes解。

34

Fisher准则的基本思路是把高维空间的点向低维空间投影，并且通过投影方向的选择，使得在(被)投影空间上，不同母体的点“尽可能分离开来”。举例来说，如有两个母体，则把m维空间的点投影到一维空间，并且通过投影方向的选择，使两个母体的投影点尽可能地分别位于直线的两侧。多个母体也是这样，首先把m维空间的点投影到r(r?m)维空间，然后选择投影方向进行分离。

当G=2时，这两种准则是等价的(正态假设下)。本章介绍Fisher准则下的两类母体的判别模型和Bayes准则下的多类判别模型；但对于Bayes准则，分别函数不同，判别函数也不同，本章主要介绍正态母体。

第二节 Fisher准则下的两类线性判别模型

一．Fisher准则的基本含义

设有A、B两个(样品)母体，若从

A中抽取nA个样品，从B中抽取nB个样品，共有

?、xp；用i=1，2，?，nA(或nB) n=nA+nB个样品；每个样品测得p个变量，即x1、x2、代表样品，j=1，2，?，p代表变量；则xAij(xBij)代表第A(或B)类中，第i个样品的第j个变量的取值。

?xA11?x?A21?...?xAnA1?x??xB11??xB21?...???xBnB1...xA1p?xA22...xA2p??.........??xAnA2...xAnAp?

xB12...xB1p??xB22...xB2p?.........??xBnB2...xBnBp??xA12两类线性判别分析的判别函数为：y?

?cj?1pjxj (2)

35

问题：如何求待定系数cj(j=1，2，?，p)？

假定判别函数已经建立，显然每个样品的p个变量代入式(2)中，就可以得到一个y值，记为：

yAi??cjxAij,yBi??cjxBijj?1j?1ppi?1,2,?,nA/nB(3)

yAi、yBi称为样品的判别计量(判别值)。 1记 yA?nA即

1y,y??AiBnBi?1nApnA?yi?1nBBi

1yA?nAyB?1nB??ci?1j?1nBpjxAijxBij(3)

??ci?1j?1j(4)称yA、yB为每类样品判别值的平均值或类平均值。

若A、B两母体存在差别，则yA与yB也会有一定的差别。使两母体分开的综合指标值y0称为两母体的分界值，或临界值(如下图所示)。

x A B

0 yA y0 yB y

显然，判别分析要求找到的判别函数y?f(x1,x2?,xp)使两类间差别愈大愈好，即

|yA?yB|或G?|yA?yB|2?max并使两类组内离差平方和愈小愈好，即

(5)

H??(yAi?yA)??(yBi?yB)2?min2i?1i?1nAnB(6)

36

结合(5)(6)得：

GI??H(yA?yB)2?(yi?1nAAi?yA)2??(yBi?yB)2i?1nB?max(7)

称(7)为Fisher准则(Fisher最大分离准则)。 二．两类线性判别函数的建立

为求极值，可使：

?I?0?cj(j?1,2,?,p)(8)

?cp，则判别函数建立。下面推导之。 从(8)解出c1、c2、对式(7)两边取对数，则LnI?LnG?LnH 故

?LnI?LnG?LnH???0 ?cj?cj?cj1?G1?H1?G?H?? 所以

G?cjH?cjI?cj?cjyA??cjxAj,yB??cjxBj

j?1j?1pp即

(9)

??G?(yA?yB)?(?cjxAj??cjxBj)?[?cj(xAj?xBj)]2

22j?1j?1j?1ppp令

dj?xAj?xBj

?cdjj?1pj则 G?()2?(c1d1?c2d2???cpdp)2(10)

?G?2(c1d1?c2d2???cpdp)dj?cj由于yAi?yA?(j?1,2,?,p)p(11)

?cj?1ppjxAij??cjxAj??cj(xAij?xAj)

j?1pj?1ppyBi?yB??cjxBij??cjxBj??cj(xBij?xBj)

j?1j?1j?1 37

H??(yAi?yA)??(yBi?yB)22i?1i?1nAnB??[?cj(xAij?xAj)]2??[?cj(xBij?xBj)]2i?1j?1i?1j?1nApnBp

∴

nbnApp?H2??2[?cj(xAij?xAj)](xAij?xAj)??2[?cj(xbij?xbj)]2(xbij?xbj)?cji?1j?1i?1j?1?2?[c1(xAi1?xA1)?c2(xAi2?xA2)???cp(xAip?xAp)](xAij?xAj)i?1nA?2?[c1(xBi1?xB1)?c2(xBi2?xB2)???cp(xBip?xBp)](xBij?xBj)i?1nB

?2c1[?(xAi1?xA1)(xAij?xAj)??(xBi1?xB1)(xBij?xBj)]??i?1i?1nAnB?2cp[?(xAip?xAp)(xAij?xAj)??(xBip?xBp)(xBij?xBj)]i?1i?1nAnB(12)令 sjj?i?1?(xAij?xAj)?2nAi?12(x?x)?BijBjnBnA?nB?2nA

为第j个变量的组内方差，

s?令 jk?(xi?1Aij?xAj)(xAik?xAk)??(xBij?xBj)(xBik?xBk)i?1nBnA?nB?2

为第j、k个变量的组内协方差(j，k=1，2，?，p)，则

?H?2(c1sj1?c2sj2???cpsjp) ?cj1?G?H?，得 I?cj?cj而

1(c1d1?c2d2???cpdp)dj?c1sj1?c2sj2???cpsjp I(j?1,2,?,p)令

1(c1d1?c2d2???cpdp)?b I(j?1,2,?,p)(13)

则 c1sj1?c2sj2???cpsjp?bdj即为：

38

?c1s11?c2s12???cps1p?bd1??c1s21?c2s22???cps2p?bd2???cs?cs???cs?bd2p2pppp?1p1(14)

式中b为常数项，它不依赖于j而变化，其大小对判别效果没有影响，故令b=1，于是

得：

c1sj1?c2sj2???cpsjp?djj?1,2,?,p解(15)式求出c1、c2、?、cp，代入(2)，即可得线性判别函数 y?c1x1?c2x2???cpxp

三．分界值计算和判别法则

当求出判别函数后，可进一步求出类平均值yA、yB，即

ny1ApA?nAi?A?yi?1?cjxAjj?1ny1Bp

B?nBi?BjB?yi?1?cjxj?1式中，xAj、xBj分别为A类和B类中第j个变量的平均值(j=1，2，?，p)，分界值y0的求法可分为以下几种情况：

（1）当母体为正态母体，且?A、?B(标准差)已知时，则 y?A0?yA?(yB?yA)?A??(16)

B（2）当nA?nB，且?A??B时，则 y0?12(yA?yB)(17)

（3）当nA与nB差别较大，?A??B时，则

yn

AyA?nByB0?n(18)

A?nBp由(18)式可导出y0?j?cjxj(19)

?1

39

(15)

式(19)中xj为两个母体中第j个变量的总平均值。 判别：

对任一个体X?(x1,x2,?,xp)，计算

Ty(x)?若y>y0，则该个体属于A类； y

j?1?cxjpj

四．判别函数的显著性检验及判别率

问题：据y判别是否有效？判别效果好坏如何？

? ? 判别函数y对区分两个 用正确判别率来衡量 母体是否显著

１．判别函数的显著性检验

1)F—检验 两个母体差异性是否显著，可用类平均yA、yB考察。 做原假设H0:yA?yB 则当H0成立时，统计量F?式中 SR?SRfR~F(1,nA?nB?2)

SDfDi?1nA?(ynAA?y)?2i?1?(ynBnBB222?y)?nA(yA?y)?nB(yB?y)

SD?i?1?(yAi?yA)?2i?1?(yBi2?yB)

fR?2?1?1,fD?nA?nB?2

给定显著水平?，F?F?(1,nA?nB?2)，否定H0，说明两组间差别显著，判别是有效的；否则无效。 注：ST?i?1?(ynAAi?y)?2i?1?(ynBBi2?y)

40

ST= SR + SD

2)马氏距离检验 通过检验母体中诸变量的类平均值之差

dj?xAj?xBj在显著差异。

设F?(j?1,2,?,p)是否足够大，即两类多元平均值在统计量上是否存

nAnB(nA?nB?1?p)D2~F(p,nA?nB?1?p)

(nA?nB)(nA?nB?2)p2式中马氏距离D?c1d1?c2d2???cpdp??cdjj?1pj(19)

在假设xAj?xBj下，F~F(p,nA?nB?1?p)

给定?，F?F?(p,nA?nB?1?p)，说明两类多元平均值间差异性显著，判别有效；否则判别无效。 ２．正确判别率

正确判别率是指属于A类的样品，根据判别值仍判别为A类母体的样品所占的百分比。若原来以A中抽取nA个样品，经过重新判别，有m个判为A类，则A类的正确判别率为

rA?n?mm?100% ?100%，错判率为eA?1?rA?AnAnA五．变量的选择

１． 根据对马氏距离的贡献挑选变量 将式(19)两边除以D得：

1?(c1d1?c2d2???cpdp)D 记

22?Dx1?c1(xA1?xB1)D2?100%?2?Dx2?cp(xAp?xBp)D?100% ????D?c(x?x)D2?100%pApBp?xpDxj(j?1,2,?,p)为第j个变量对马氏距离的贡献，挑选贡献最大的变量参加建立判别函数。

２． I值最大法挑选变量

Fisher最大分离准则

GI??H(yA?yB)2[?(yAi?yA)2??(yBi?yB)2]i?1i?1nAnB?max

41

若考虑一个变量xj，则由上可得每个变量的Ij值(j=1，2，?，p)为：

Ij?(xAj?xBj)2[?(xAij?xAj)2??(xBij?xBj)2]i?1i?1nAnB

显然Ij越大的变量，分辨母体的能力越强。

注：只考虑变量独立影响，当变量不独立时，用逐步判别分析。

第三节 Bayes准则下的多类线性判别模型

一、几个概念

1) 先验概率 qg (g=1,2,…,G)

未经观测之前，根据已有资料知道某一个体来自母体g的概率。因为某一个体只能属于其中一类(互不相容)，故根据概率加法定理

g?1?qGg?q1?q2?...?qG?1

2) 条件概率p(x|Ag) (g=1,2,…,G)

在已知个体x来自母体Ag的条件下观测到的个体x的概率。当分别函数为连续函数时，可看作Ag的概率密度函数。

3) 后验概率p(Ag | x) (g=1,2,…,G)

在已知观测到个体x的条件下，个体x来自母体Ag的概率(个体x属于的Ag概率)。亦即在个体x出现的情况下，x来自Ag的概率。 先验概率和后验概率的区别在于：

前者是指对于个体未经观测之前知道的该个体来自母体Ag的概率。 后者是指已经对该个体进行观测之后知道的该个体来自母体Ag的概率。

根据上述定义，Bayes公式可写成

p(g|x)?qgp(x|Ag)i?1?qiG?qgp(x|Ag)p(x)

ip(x|Ai)(x)?后一步由全概率公式pi?1?qGp(x|Ai)得到。

(g|x)?分布函数为连续函数时：pqgfg(x)Giii?1

?qf(x) 4) 错分概率p(h | g)

42

个体实属Ag但把它划为Ah的概率。 任何一种划分都存在着错分现象。

如果我们把个体X?(x1,x2,?,xp)看成是p维欧氏空间R上的一个点，Bayes准则设法将R划分成G个互斥的子空间R1，R2 ，… ,Rg，如果样品落到Rg，就把该个体归为Ag 。

故 p(h|g)?RhTdx ?f(x)g

5) 错分损失L(h|g)

把原属于Ag的样本划归Ah母体中的损失。

错分损失是一个具体的量，但是在实际工作中很难估计。 约定 L(h|g)=0 g=h L(h|g)>0 g?h

进一步约定L(h|g)=L(g|h)(=1) g?h

二、Bayes准则

若已知任一个体来自A1，A2 ，… ,Ag 共G个母体，Bayes判别法将样本空间R划分成R1，R2 ，… ,Rg共G个互不相交斥的子空间。对于R的任一分法，都存在着错分现象。假设fg (x)和qg已知，则将Ag错分的平均损失为

h?1h?g(h|g)p(h|g) g=1,2,… ,G ?Lp(h|g)?dx ?f(x)gGRh对应于一个分法的总平均损失

WR＝?qg?L(h|g)p(h|g)

g?1h?1h?gGG使上述损失最小的分法就称为Bayes解。

对于具体样品，记Eh (x)为把x划入h类所造成的损失，则

Eh(x)?g?1g?h?GL(h|g)p(h|x)?qg?Gqgfg(x)i?1g?1g?h?qf(x)iiGL(h|g)

43

可以证明，若Eh*(x)?min{E(x)}

h1?h?G则把x划归Ah* 类，这样的分法就是一个Bayes解。 Bayes准则可归纳如下：

假设各母体Ag (g=1,2,…,G)的概率密度fg (x)已知，先验概率qg (x)已知，并且错分损失相等，这时可建立判别函数

qg fg (x) g=1,2,…,G 如果 qg*(x)fg*(x)?f(x)} max{q(x)gg1?g?G则把个体x划归母体g*。

该准则对任意分布都适用。

三、正态母体判别函数的建立和样品判别方法 1. 线性判别函数

一元正态分布N(?,?2)的概率密度函数为

21(x??)f(x)?exp[?] 22??2? G个母体Ag(g=1,2,…,G)都服从多元正态分布N(ag,?)(假设?与g无关)，其概率密度函数

为

fg(x)?Σ1?12p2exp[?(2?)1T?1(x?ag)Σ(x?ag)] 2 其中 X?[x1,x2,...,xp] 是一个p维矢量； ag?[ag1,ag2,...,agp] 是数学期望矢量；

TTΣ＝[?kj]p?p??11??21???...???p1?12...?1p??22...?2p??

.........???p2...?pp?T为协方差矩阵；?-1是?的逆矩阵，|?-1|是?-1的行列式。 在实际问题中，用xg?[xg1,xg2,...,xgp]

xgj作为ag的估计值；用S＝[skj]p×p

1?ngi?1?xnggij

skj?g?1i?1??(x?1Gnggij?xgj)(xgik?xgik)/(n?G)

作为?的估计值。记S的逆矩阵S

?[skj]p?p，则此时正态母体的密度函数可改写为

44

fg(x)?S1?12p2exp[?(2?)1T?1(x?xg)S(x?xg)] 2判别函数为qgfg(x)。由于我们只关心寻求使qgfg(x)最大的g,故可取对数

??11??S2?1T?11T?1T?1ln(qgfg(x))?ln(qg)?ln??xSx?xSx?xgSxg ?gp22?(2?)2???去掉与g无关的项，并记

T?1Cg?xgs?[cg1,cg2,...,cgp]

cgj?k?1?spkjxgk j=1,2,…,p; g=1,2,…,G

cg0最后可得判别函数

11??xgTS?1xg??22j?1?cpgjxgj

yg(x)?lnqg?cg0?Tj?1?cpgpxj g=1,2,…,G

把个体X?[x1,x2,...,xp]的值代入判别函数，计算出y1(x), y2(x),…, yG(x), 如果 yg*(x)?则把x划归第g*个母体。

2. 计算步骤 1) 读入数据

max{y(x)}

g1?g?G?x111x112?xx122?121?......??x1n11x1n12?x211x212??x221x222x??......??x2n21x2n22?......??xG11xG12??xG21xG22?......??xGnG1xGnG2?

45

...x11p?...x12p??......??...x1n1p?...x21p??...x22p?......?

?...x2n2p?......??...xG1p??...xG2p?......??...xGnGp??

2) 求平均值

xgj1?ngi?1?xnggij

xg?[xg1,xg2,...,xgp]T

3) 求协方差矩阵的估计S及S-1 skj?g?1i?1??(xGnggij?xgj)(xgik?xgik)/(n?G)

n?g?1?nGg

S?[skj]p?p S?1?[skj]p?p

4) 求判别函数系数

T?1Cg?xgs?[cg1,cg2,...,cgp]

cgj?k?1?sppkjxgk

cg01??2j?1?cgjxgj

判别函数 yg(x)?lnqg?cg0?j?1?cpgpxj

qg?ng以样本频率代替?n ??1?? 各类相等

G??qg 人为给定??g无论哪种，都应有

g?1?qG?1

5) 对样品进行判别(计算后验概率)

(Ag|x)?公式 peGh?1yg(x)

yh(x)?e 46

对每一个样品，先计算yg(x) g=1,2,…,G

挑出最大的g*,则判别样品属于g*类，后验概率

p(Ag*|x)?eGh?1yg*(x)?yh(x)1h?1?e?eG

y?h(x)(x)?yh(x)?yg*(x) (可防止溢出) 其中 y?h

四、判别效果检验(略)

第四节 逐步判别分析方法

一．逐步判别分析的基本思想

逐步判别分析和逐步回归分析的基本思想相似，它们都是根据每一个变量在各类(组)判别式中所起的判别作用的重要性不同来挑选判别效果最好的变量进入判别式，同时从判别式中剔除那些由于新变量的引入而失去判别作用的变量，使最后的判别式中，只保留对母体判别能力较强的变量。

逐步判别分析的做法是：

（1）比较变量判别能力，挑选变量，检验显著性，若显著则引入。 （2）再把未选入的变量同已选入的变量结合比较判别能力，从中选择分辨能力最大的变

47

量，检验显著性，若显著则引入。 （3）当第二个变量选入后，考察第一次引入的那个变量是否由于新变量的引入而判别能

力下降。故检验第一个变量在第二个变量存在时的判别能力是否显著。如果不显著，则剔除；否则，再考虑引进另一个变量。再把未选中的每一个变量与已选中的两个变量组合，重复上述步骤，直至既无变量引入，又无变量剔除为止。 （4）最后，利用选入的变量建立判别函数。 二．变量的判别能力及变量取舍的标准 １．变量的综合判别能力

假设收集了G类样品，第g类有ng个样品(g?1,2,?,G,共有n?，?n个样品）gg?1G其组内离差平方和矩阵和总离差平方和矩阵为一阶矩阵，只有一个元素，即：

W1?w11???(xgj?xg)2

g?1i?1ngGngT1?t11???(xgi?x)2

g?1i?1G则 U1?|W1|表示变量x的判别能力(U为Wilks统计量)。 |T1|如有两个变量(如x1、x2)时，

?wW2??11?w21其中

w12?w22???tT2??11?t21t12? t22??Gng??wkj???(xgik?xgk)(xgij?xgj)g?1i?1??Gng?t?(xgik?xk)(xgij?xj)?kj??g?1i?1?(k,j?1,2)

则 U1,2?|W2||T2|(1)

反映了两个变量x1,x2组合在一起分辨母体的能力。即U1,2值越小，则x1,x2的判别能力越强。

当有L个变量(如x1,x2,?,xL)时，其组内变差矩阵与总变差矩阵为L阶矩阵，即：

?w11?wWL??21????wL1

w12w22?wL2?w1L??w2L???????wLL??t11t12?ttTL??2122?????tL1tL248

?t1L??t2L?? ?????tLL?其中

Gng??wkj???(xgik?xgk)(xgij?xgj)g?1i?1??Gng?t?(xgik?xk)(xgij?xj)?kj??g?1i?1?(k,j?1,2?,L)

则 U1,2,?,L?|WL||TL|(2)

反映了L个变量组合在一起分辨母体的能力。即U1,2,?,L值越小，则变量的判别能力越强。 ２．未选变量的判别能力及引进变量的标准

1）未选变量的判别能力 假定变量x1,x2,?,xL已经给定，然后再添加一个新变量

xr(r?L)，现在讨论xr的判别能力

把L+1个变量分为两组：第一组是前L个已给定的变量；第二组仅有xr。那么L+1个变量的组内变差矩阵为

WL?1?w11?w?21?????wL1?w(L?1)1???w12w22?wL2????w1Lw2L?wLLw(L?1)2?w(L?1)Lw1(L?1)?w2(L?1)?????

wL(L?1)?w(L?!)(L?1)????W12??WWL?1??11，其中 ??W21W22?W11?WLW12?W21?[w1(L?1),?,wL(L?1)] Wrr?W(L?1)(L?1)?wrr将行列式进行变换，可得：

|WL?1|?|W11|1?1W11W12?1Wrr?W21W11W12

?1?W11Wrr?W2W11W12(l)令 wrr?Wrr?W21W11W12，则|WL?1|?|W11|wrr

(l)(l)?1同理 |TL?1|?|T11|trr

49

∴

WL?1TL?1(l)WLwrr? (l)TLtrr故

U1,2,?L,r?WL?1TL?1?U1,2,?LUr|1,2,?,L

其中U1,2,?L,r表示变量x1,x2,?,xL及xr的判别能力；

U1,2,?L表示前L个变量x1,x2,?,xL的综合判别能力；

Ur|1,2,?L表示在给定前L个变量x1,x2,?,xL的条件下，变量xr的判别能力，它表示xr与x1,x2,?,xL结合在一起时的判别能力，记为Ur|(L)。

2）引进变量的标准 假设已经计算了l步，并引进了L个变量(包括l=0；L=0)，尚有

xL?1,xL?2,?,xL?m个变量尚未引入，现在要确立第l+1步及再引入一个新变量

xr(r=L+1，L+2，?，L+m)的标准。

在给定x1,x2,?,xL的条件下，每一个未选入的变量xr的判别能力均可由

Ur|(L)(l)wrr?(l)trr(r?L?1,L?2,?,L?m)(3) 给出。

其中，必有一个最小者，记为 Ur|(L?minUr|(L)**??(r?L?1,L?2,?,L?m)

*说明变量xr的判别能力最强，用与Ur|(L)等价的F近似式

1?Ur*|(L)n?G?LF(G?1,n?G?L)?~F(G?1,n?G?L)*G?1Ur|(L)*(4)

进行显著性检验。若F?F?，则认为xr判别能力显著，应该引入判别函数。(4)式中F作为引入变量的标准，称为“引入F”。

３．已选变量的判别能力及剔除变量的标准

1）已选变量的判别能力 在逐步引入变量时，应重新估计已选入的变量。假设已经计算了

l步，并引入了包括xr在内的L个变量，现要确定第l+1步剔除变量xr的标准。记 Ur|(L?1)(l?1)wrr?(l?1)trrr?1,2,?,L(5)

即前l-1步引入了不包括xr在内的L-1个变量，Ur|(L?1)值愈大，xr的判别能力愈弱。

50

本文来源：https://www.bwwdw.com/article/1r0o.html