数学地质实验指导书（教材）

《数学地质》

实 验 指 导 书

二零零六年八月

说 明

一、该指导书所属课程：

《数学地质》 二、适用专业：

地质工程 三、实验总学时数：

20学时。可根据需要增开适当的课外机时。 四、各实验项目时数：

实验一 预处理与统计 2学时实验二 线性回归分析 4学时实验三 多元线性回归分析 4学时实验四 趋势面分析 4学时实验五 聚类分析 4学时实验六 两类判别分析 4学时实验七 贝叶斯多类判别分析 4学时实验八 有序地质量最优分割法 2学时五、前修课程知识 1． 计算机编程语言； 2． 线性代数； 3．

概率统计。

必做实验选做实验必做实验必做实验必做实验必做实验选做实验必做实验 目 录

实验一 预处理与统计 ................................................................................................................... 1 实验二 一元线性回归分析 ............................................................................................................ 7 实验三 多元线性回归分析 ...........................................................................................................11 试验四 趋势面分析 .................................................................................................................... 15 试验五 聚类分析 ........................................................................................................................ 20 试验六 两类判别分析 ................................................................................................................. 24 试验七 贝叶斯多类判别分析 ...................................................................................................... 28 试验八 有序地质量最优分割法 ................................................................................................... 32

实验一 预处理与统计

一、目的：

通过完成数据统计和预处理程序的设计和实现及完成算例，掌握统计一组数据的极值、均值、方差、变异系数及进行数据预处理的方法。

二、方法概要

1、进行统计和预处理的原因、目的和应注意的问题

（1）原因

原始数据可能有强非对称性，存在孤立值，大多数的统计方法应用原始数据时存在大而且不是偶然的残差等问题，通过改变表达方式，有时可以增强信息的显示，而这种改变不仅需要改变数值的单位，而且可能改变数据的基本测量尺度；

（2）目的

? 使变量尽可能为正态分布（如回归分析要求因变量为正态分布，要求自变量和因变量之间

具有足够的相关关系）； ? 统一变量的数据尺度；

? 使变量之间的非线性关系转换为线性关系；

? 用新的数目少的相互独立的变量代替相互联系的原始变量； ? 方便用简单自然的方式进行解释； ? 帮助理解数据的特征。 （3）注意问题

? 数据范围：只有数据变化范围相对较大，变换才显著；

? 变换是很重要的工作，变换不当则适得其反；所以在认真研究分析的基础上进行，有时要

通过多次试验才能找到合适的变换方法；

? 有些行业中，有些强制性变换或习惯使用的变换，工作中应遵循；

? 变换后数据的可解释性也很重要，有时为了不影响解释，宁可不对其转换。

2、鉴别并剔除异常值

（1）3-σ法则（拉依达准则） 由于随机误差服从正态分布规律，因此

P[|?|?3?]?P[x???3?]?99.7% （1-1）

由式1－1知，误差|?|?3?出现的概率只有0.3%,也就是说,在1000次测量中,误差大于3?的情况只可能出现3次。因此，在有限次测量中，若某次测量值的误差大于3σ时，则认为该测量值

1

含有过失误差，应予以舍去，这就是3?法则。

一般的测量时误差?是得不到的，只能得到残差v?x?x，而总体标准差?也是得不到的，故?，??通常取为样本标准差，只能用它的估计值 ?即??=S。因此，3σ法则只能按下述规则实际应用：

对于测量数据x1,x2,?,xn，若某个测量值xi（1?i?n）对应的残差满足

vi?xi?x?3S （1-2）

则将xi舍去。

根据3σ法则对实验数据进行处理，也有犯“弃真”的错误。就是将一些误差较大但并不含有过失误差的正常测量值当作含有过失误差的测量值舍去了。3σ法则“弃真”概率很小，且随着测量次数的增加而减少，最后稳定与0.3%。一般要求n>10。当n≤10时，即使测量数据中含有粗大误差，用3σ法则也不能判别出来。

（2）Dixon准则

第一步 将样本从小到大顺序排列，得次序统计量：

x(1),x(2),?,x(n)

xmin?x(1),xmax?x(n)称为极端值。

第二步 对不同的n求极端值，选择计算不同的统计量。

Dixon[x]?x(n)?x(n?1)x(n)?x(2)(n) （1-3）

第三步 对比，若计算出的统计量>临界值，则认为相应的极端值为异常值。

3、统计

（1）平均值

设有一批样本数据（x1 ，x2 ，?，xn），其平均值为x

x?1n?xi （1-1）

ni?1（2）样本方差

设有一批样本数据（x1 ，x2 ，?，xn），其离差平方和的平均值为样本方差，记S2

S2?1n?1i?1?(xi?x) （1-2）

n2方差反映了数据的离散程度，其值越大数据越分散；其值越小数据就较多的集中在平均值附近。但它是有量纲的，受到量纲、量级的制约。

（3）变异系数及正态分布检验

样本数据的标准差于平均值之比为变异系数，记Cr

2

Cr?Sx100% （1-3）

变异系数它是无量纲的，克服了量纲的影响，给出数据相对变化性度量，能较好的反映出数据变化程度的大小。

当变异系数大于0.5，一般是非正态总体；正态总体的变异系数小于0.33；当变异系数介于0.33～0.5之间时，总体多数服从对数正态分布。

4、统一量纲

由于测量数据单位不同，或数量级不同，多种变量观测值的变化范围往往很大，如果用原始数据计算，可能突出某些数量级特别大的变量作用，压低了数量级小的变量的作用，使计算失去意义。为使变量有相等的权，必须对原始数据进行变换，使其达到量纲统一。

（1）标准化变换

x?'ijxij?xjSj?xij?xj1n?1i?1?(xij?xj)n2 （1-4）

式中xij为原始观察值，xj为第j个变量的算术平均值，Sj为第j个变量的标准差。i=1，2，?，n为样品号，j=1，2，?，m为变量号。

变换后，变量平均值为0，方差为1，各变量处于同一量纲，每个变量在变化前后的相关程度不变。按照几何意义，标准变换相当于将坐标原点移至中心（平均数）位置。这种变换适合于量纲不同和数量大小不均一的连续型原始数据，如化学分析数据。

（2）极差变换

x'ij?xij?xjminxjmax?xjmin （1-5）

式中xij为原始数据，xjmin为第j个变量的最小值，xjmax为第j个变量的最大值。i=1，2，?，n为样品号，j=1，2，?，m为变量号。

变换后，数据处于同一量纲，其最大值为1，最小值为0，所有数据在0~1之间变化。变化前后的相关程度不变。按照几何意义，相当于把坐标原点移至变量最小值的位置上。极差变换适合于量纲和大小不一样的连续型原始数据的变换。

（3）均匀化变换

xij?'xijxj （1-6）

式中xij为原始数据，xj为第j个变量的算术平均值。变换后数据量纲一致，数据都是在1附近波动的相对数，数学期望为1，变量与平均数之差的期望为0。此变换适合于比例型数据，如长度、

3

体积、质量等。

5、正态变换

正态变换是指把呈非正态分布的数据变换成正态或近似正态分布的数据变换。 （1）角变换

角变换的作用是使左偏和右偏的不对称分布变成近似正态分布，方法是把原始数据变为0°～90°之间的角度。公式：

?'x?arcsin??i??x'?arccosi??xi10 （1-7）

xi10nn式中xi'为变换后的数据；xi为原始数据；n=1，2，?为正整数，取变量原始观察值的最大值的整数位数。变换前后和其它变量的相关性略有差异。

（2）平方根变换

使右偏不对称分布变为正态分布。公式：

xi?'xi?C （1-8）

C为常数。

平方根变换适用于服从泊松分布的离散型变量。如露头个数，距主断裂带的距离等。 （3）对数变换 公式：

xi?log(xi?C) （1-9）

'C为常数。

对数变换适用于服从对数正态分布的数据。如化探数据及有色、稀有、重金属的品位数据等。

三、程序设计框图

预处理和统计模块的流程图如图1-1所示。

后续的不同多元统计方法对原始数据的分布有不同的要求，但对原始数据都需要进行预处理和统计，所以建议统一输入文件的格式，预处理和统计模块处理后输出文件的格式也按该定义的格式定义。

文件的格式建议为：

样品个数，变量个数

第1个样品号，该样品第1个变量的观测值，第2个变量的观测值，… 第2个样品号，该样品第1个变量的观测值，第2个变量的观测值，… 第3个样品号，该样品第1个变量的观测值，第2个变量的观测值，… …

不是所有的原始数据都进行系统的统计和预处理，程序应该提供强大和灵活的人机交互界面，

4

方便用户根据实际情况和需要，选择不同的数据变换方法。

开始打开标准格式的原始数据文件读入样品个数n,变量个数m，原始数据矩阵x(i,j)形成新的样品个数n?,变量个数m?，原始数据矩阵x?(i,j)数据统计模块求平均值数据统一量纲模块求平均值求样品方差计算平均值x(j)?1n?x?x,j?ni?11n计算样品方差：S(j)??(x(i,j)?x(j))2和S(j)n?1i?13-σ法则2求样品方差求变异系数和正态分布检验i=1,j=1计算残差v(i,j)?x(i,j)?x(j)是否正态标准化变换极差变换均匀化变换选择正态变换的方法YV(I,j)>3SY删除样品I和变量JN形成新的样品个数n??,变量个数m??，原始数据矩阵x??(i,j)角变换对数变换平方根变换j=j+1j≤mNJ=1i=i+1输出变换后的标准格式数据文件开始Yi≤n

图1-1 数据预处理和统计模块流程图

四、算例

（1）请对以下15个样品数据进行数据统计和数据预处理

样品 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 1.766115 1.728760 1.715418 1.459392 1.738305 1.713323 1.751741 1.853637 1.649627 1.720903 1.737908 1.631849 1.720407 1.828789 1.709948 2 3.481379 3.234192 3.682390 3.982461 3.505710 3.293934 3.482815 3.461214 3.773592 3.392639 3.471311 3.824918 3.292416 3.497142 3.675595 3 52.62 52.52 61.10 30.36 71.90 53.58 66.16 64.26 55.80 51.70 59.20 56.12 56.06 72.88 67.90 变量 4 5 5.70 15.27 2.20 38.92 5.20 23.35 33.59 28.21 4.00 14.59 1.80 37.84 4.00 18.59 5.90 14.27 7.80 23.89 2.60 37.19 3.70 30.48 10.50 17.3 4.60 36.00 6.05 14.49 5.40 15.84 6 9.95 2.40 2.78 3.09 2.33 3.23 4.46 5.56 4.19 3.36 2.75 8.37 1.17 3.09 4.39 7 5.82 1.04 2.86 0.89 1.04 1.09 1.38 2.27 3.11 2.44 1.18 2.07 1.23 1.04 2.96 8 9.60 1.40 1.36 1.58 1.55 1.80 1.61 3.98 2.23 1.00 0.94 4.90 0.31 1.28 2.19 9 1158 1578 1372 1383 1285 1564 1320 1230 1350 1549 1479 1233 1555 1290 1282 5

（2）请对以下15个样品数据进行数据统计和数据预处理 样品 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 1.766115 1.72876 1.715418 1.459392 1.738305 1.713323 1.751741 1.853637 1.649627 1.720903 1.737908 1.631849 1.720407 1.828789 1.709948 2 3.481379 3.234192 3.68239 3.982461 3.50571 3.293934 3.482815 3.461214 3.773592 3.392639 3.471311 3.824918 3.292416 3.497142 3.675595 3 52.62 52.52 61.1 30.36 71.9 53.58 66.16 64.26 55.8 51.7 59.2 56.12 56.06 72.88 67.9 4 5.7 2.2 5.2 33.59 4 1.8 4 5.9 7.8 2.6 3.7 10.5 4.6 6.05 5.4 变量 5 15.27 38.92 23.35 28.21 14.59 37.84 18.59 14.27 23.89 37.19 30.48 17.3 36 14.49 15.84 6 9.95 2.4 2.78 3.09 2.33 3.23 4.46 5.56 4.19 3.36 2.75 8.37 1.17 3.09 4.39 7 5.82 1.04 2.86 0.89 1.04 1.09 1.38 2.27 3.11 2.44 1.18 2.07 1.23 1.04 2.96 8 9.6 1.4 1.36 1.58 1.55 1.8 1.61 3.98 2.23 1 0.94 4.9 0.31 1.28 2.19 9 1158 1578 1372 1383 1285 1564 1320 1230 1350 1549 1479 1233 1555 1290 1282

6

实验二 一元线性回归分析

一、目的：

通过对一元线性回归分析程序设计及完成算例，掌握一元线性回归分析的基本原理和方法。

二、方法概要

1、原始数据预处理

对n?对原始数据（xi,yi）(i=0，1，?，n?)数据进行预处理（见实验一）。

2、数学模型建立

（1）作散点图

一元线性回归只能处理两个变量之间的关系，它又被称为直线拟合，假设x为自变量，y为因变量，对预处理后n对数据（xi,yi）(i=0，1，?，n；其中n?n?)如果把各个点数据画在坐标纸上（作散点图），各点近似分布呈一条直线，则可用一元线性回归。

（2）选择数学模型

??a?bx （2-1） yy

z

散点图

（3）参数a, b的估计

?为根据回归方程的得到的因变量y的计算值；a，b为回归方程中的系数；x为自变量。由于y测定结果中不可避免会有实验误差，因此用最小二乘法的原理估计回归直线中的系数a，b。

假设实验得到了n对数据（χ

i ，

??a?bx，则对于xy）(i=0，1，?，n)并得到了回归方程y的一系列变量x1 ，x2 ，?，xn，根据回归得到方程可得到因变量的一系列计算值

?1?a?bx1y?2?a?bx2y???n?a?bxny

?1(i=0，1，?，n)之间都有偏差，每一个实测值yi(i=0，1，?，n)与它相对应的一个计算值y也可称残差，计算公式：

?i?yi?(a?bxi) （2-2） vi?yi?y所有测试数据的残差平方和为

7

nn2iniQe??vi?1??[yi?1?]??y2?[yi?1i?(a?bxi)]2 （2-3）

如果回归方程y??a?bx是合理的a，b为最佳。则所得到的残差平方和应达到最小值。 得

n??xiyi?nxyLxy??1?b??in?22 （2-4） Lxx?xi?n(x)?i?1?a?y?bx?列出回归方程

y?a?bx （2-5）

3、回归方程显著性检验

（1）将自变量n组试验数据依次代入回归方程（式3－8），求出因变量n个回归值

?i?a?bxi （i?1,2,?,n） （2-6） y（2）计算S总，S回和S剩：

nS总??yi?1n2i2?ny （2-7）

S回???yi?12i?ny （2-8）

2S剩＝S总?S回 （2-9）

（3）计算F统计量

F?S回1S剩n?2 （2-10）

（4）计算相关系数r

r?LxyLxx?Lyy （2-11）

（5）查表检验

给定显著水平α，查表得F?(p,n?2),r?(n?2)。若F?F?(p,2)或r?r?(n?2)，则证明回归方程显著；否则，回归方程无实用价值。

4、对因变量进行区间估计

（1）估计剩余标准值

???S剩(n?2 （2-12）

8

（2）求出因变量估计值

当自变量取定值x时，由回归方程可求得 因变量数学期望的估计值y?。 （3）写出预测区间

??2??,y??2??） 95%置信区间:（y??3??,y??3??） 99.7%置信区间：（y三、程序设计框图

建议对原始数据以以下格式输入：

样品个数，2

第1个样品号，该样品第1个变量的观测值，第2个变量的观测值 第2个样品号，该样品第1个变量的观测值，第2个变量的观测值 第3个样品号，该样品第1个变量的观测值，第2个变量的观测值 …

首先要对数据进行预处理，处理完成后，形成以上格式新的文件（详见图1－1）。以下步骤如流程图（图2－1）。

开始打开标准格式的原始数据文件数据预处理和统计模块输出变换后的标准格式数据文件F>Fα且r>rαY输出F,r,S，F?,r???a?bx读入待预测样品x值，求出y读入n对数据（xi，yi）（I＝1，2，?，n）计算x,y,lxx,lxy,lyy(见式2－3）计算a,b(见式2－3）计算剩余标准差???S剩/(n?2)?i?a?bxi(i?0,1,?，n)求出yN??2??L21＝y??2??L22＝yn?yi2?ny2?S总??i?1n??i2?ny2?S回?y计算?i?1?S?S?S剩总回?S回1?F???S剩n?2???3??L31＝y??3??L32＝y??计算r?LxyLxx?Lyy输出95％的置信区间（L22,L21)和99.7%的置信区间（L32,L31)结束输入F?,r?（或者建立数据库，输 入?，系统自动查询F?,r? ）

图2-1 一元线性回归分析模块流程图

四、算例

（一）从某煤矿采集11个煤样，分别测定煤的发热量（QDT）和煤灰分（Ag）含量，获得如下表数据：

g 9

样品号 g1 2 3 7.7 4 7.2 5 6.8 6 6.2 7 5.6 8 9 10 11 QDT（千卡） 8.30 7.8 g5.5 5.00 4.70 4.3 0.03 0.05 0.08 0.10 0.15 0.20 0.27 0.30 0.34 0.40 0.45 A g试建立煤的发热量（QDT）对煤灰分（Ag）的回归方程，并检验该回归方程的显著性。

（二）煤矿脉中13个相邻样本点处某种伴生金属的含量数据如下表：

样品号 距离x 含量y 样品号 距离x 含量y 1 2 106.42 8 11 110.59 2 3 108.20 9 14 110.60 3 4 109.58 10 15 110.90 4 5 109.50 11 16 110.76 5 7 110.00 12 18 111.00 6 8 109.93 13 19 111.20 7 10 110.49 试建立y对x的回归方程，并进行显著性检验。 （三）收集到广西晚二迭世含煤建造30个工程点煤系厚度（H）和煤层厚度（M）的资料，如表三（提示：用抛物线方程y?a?bx?cx2）作二元一次线性回归。

样点编号 1 煤厚（M） 1.80 煤系厚（H） 121.2 样点编号 11 煤厚（M） 0.40 煤系厚（H） 528.5 样点编号 21 煤厚（M） 3.35 煤系厚（H） 240.9 2 1.05 57.3 12 1.10 480.0 22 2.06 154.8 3 0.56 63.1 13 3.49 420.0 23 1.63 116.7 4 0.0 38.0 14 0.46 418.5 24 3.70 183.7 5 2.30 154.6 15 0.0 20.0 25 1.18 127.0 6 0.92 138.8 16 2.80 99.1 26 1.8 165.4 7 2.50 134.6 17 2.54 378.5 27 4.31 192.4 8 2.93 126.9 18 0.92 596.3 28 2.22 179.1 9 3.26 196.1 19 4.33 155.8 29 2.70 230.0 10 2.28 328.7 20 2.84 357.7 30 4.00 226.0 试建立煤层厚度（M）对含煤建造厚度（H）的回归方程，并作显著性检验。

10

实验三 多元线性回归分析

一、目的：

通过对多元线性回归分析程序设计及完成算例，掌握多元线性回归分析的基本原理和方法。

二、方法概要

设有自变量x1,x2,?,xp，因变量y，共做n次实验。若y与x1,x2,?,xp间有线行关系，回归方程则为：

y?b0?b1x1?b2x2???bpx （3-1）

显而易见，只要确定了b0,b1,?,bp各回归系数，方程也就确定了。

1、建立原始数据矩阵，并进行必要的预处理

（1）确定因变量与自变量，形成n行p?1列的矩阵：

?x11?x21??x??n1x12x22xn2????x1py1?x2py2?? （3-2） ?xpyn??（2）进行预处理（详见图1－1）

2、数学模型建立

设经过预处理的原始数据，表示如下

?x11x12?x1p??y1??xx?x2p??y? Y?2 （3-3） X??2122??????xx?x???yn??p??n1n2?的偏差平方和Q达到最小，即由最小二乘法知道，b0,b1,?,bp使得全部观察值y与回归值y使

nniQ??(yi?1?i)??y2?(yi?1i?b0?b1xi1?b2xi2???bnxip)?最小 （3-4）

2根据微积分学中的极值原理，b0，b1，?，bn，应是下列方程组：

n???i)?0Q??2?(yi?y???b0i?1 (i?1,2,?,n；j?1,2,?,p) （3-5） n???i)Xij?0?Q??2?(yi?y?i?1??bj3－5式是求解回归系数的正规方程组，可以写成以下形式

11

?XX??B?B??XX??X?Y

?1 （3-6）

X?Y （3-7）

求解正规方程组，得到各回归系数b0,b1,?,bp。 列出回归方程

y?b0?b1x1?b2x2???bpxp （3-8）

3、回归方程显著性检验

（1）将自变量n组试验数据依次代入回归方程（式3－8），求出因变量n个回归值

?i?b0?b1xi1?b2xi2???bpxip （i?1,2,?,n） （3-9） y（2）计算S总，S回和S剩：

nS总??yi?1n2i2?ny （3-10）

S回???yi?12i2?ny （3-11）

S剩＝S总?S回 （3-12）

（3）计算F统计量

F?S回pS剩n?p?1 （3-13）

（4）计算相关系数r

r?S回S总 （3-14）

（5）查表检验

给定显著水平α，查表得F?(p,n?p?1),r?(n?2)。若F?F?(p,n?p?1)或r?r?(n?2)，则证明回归方程显著；否则，回归方程无实用价值。

4、对因变量进行区间估计

（1）估计剩余标准值

???S剩(n?p?1) （3-15）

（2）求出因变量估计值

?。 当自变量取定值x1,x2,?,xp时，由回归方程可求得因变量数学期望的估计值y（3）写出预测区间

??2??,y??2??） 95%置信区间:（y??3??,y??3??） 99.7%置信区间：（y 12

三、程序设计框图

建议以实验一的原始数据格式输入，首先要对数据进行预处理，处理完成后，形成实验一的原始数据格式的新文件（详见图1－1）。以下步骤如流程图（图3－1）。

开始打开标准格式的原始数据文件数据预处理和统计模块输出变换后的标准格式数据文件读入n对数据（xij,yi）（i＝1，2，?，n;j＝1，2，?，p）,形成原始矩阵（见式3－3）?1计算X?,(XX?)?1,?XX??X?Y(见式3?7）请参考线性代数矩阵?1根据B?XX??X?Y,求得b0,b1,?,bp运算。?i?b0?b1xi1?b2xi2???bpxip求y（i?1,2,?,n）n?S?yi2?ny2?总?i?1n??i2?ny2?S回?y计算?i?1?S?S?S剩总回?S回p?F???S剩n?p?1???计算r?S回S总输入F?,r?（或者建立数据库，输 入?，系统自动查询F?,r? ）F>Fα且r>rαY输出F,r,S，F?,r???b0?b1x1?b2x2???bpxp读入待预测样品x1,x2,?,xp值，求出y计算剩余标准差???S剩/(n?p?1)N??2??L21＝y??2??L22＝y??3??L31＝y??3??L32＝y输出95％的置信区间（L22,L21)和99.7%的置信区间（L32,L31)结束

图3-1 多元线性回归分析模块流程图

13

四、算例

（一）、某矿区从18个矿样中测得Cu,Pd,Pt的含量如下所示，试求Pt对Cu和Pd的回归方程，并检验其有无实用价值。 编号 1 2 3 4 5 6 7 8 9 x1 Cu(%) 0.36 0.05 0.24 0.05 0.22 0.13 0.10 0.07 0.16 x2 Pd(g/t) 0.169 0.035 0.248 0.015 0.216 0.086 1.300 0.009 0.160 x3 Pt(g/t) 0.176 0.020 0.262 0.035 0.234 0.092 1.669 0.006 0.150 编号 10 11 12 13 14 15 16 17 18 x1 Cu(%) 0.06 0.53 0.06 0.21 0.03 0.33 0.01 0.21 0.08 x2 Pd(g/t) 0.023 0.329 0.065 0.283 0.030 0.21 0.030 0.252 0.053 x3 Pt(g/t) 0.045 0.292 0.080 0.235 0.030 0.237 0.050 0.24 0.103 （二）、由于碳、氢、氧是煤燃烧过程中产生热量的主要元素，对某矿（褐煤、常烟煤、肥煤、

r焦煤和无烟煤）中取12块煤样，经分析化验后，其发热量QDT（焦耳/克）与Cr(%)、Hr（%）、Or（%）

元素的数据如下。问：今后能否不再对该矿煤进行发热量测定，而由元素分析结果对其进行预测？

编号 1 2 3 4 5 6 x1 62 70 75 75 78 80 x2 5.0 6.0 6.5 5.0 5.5 5.2 x3 15 20 25 10 15 4.0 y rDT rrC(%) H（%） Or（%） Q（J/g） 6.5 7.0 7.2 7.5 7.7 8.0 编号 7 8 9 10 11 12 rrr C(%) H（%） Or（%） QDT（J/g）x1 x2 x3 y 85 88 90 90 92 95 6.0 3.0 5.0 2.5 3.0 3.5 6.0 3.0 5.0 2.5 3.0 3.5 8.4 8.5 8.8 8.0 8.5 8.7

14

试验四 趋势面分析

一、目的：

通过趋势面分析程序设计及完成算例，掌握趋势面分析方法原理及工作步骤。

二、方法概要

1、建立原始数据矩阵，并进行必要的预处理

（1）选取应变量

选取什么样的变量进行趋势分析，取决于研究对象和研究目的。 （2）确定控制点

在确定控制点时，一般要考虑以下几个方面：①控制点的代表性、真实性和可靠性。②控制点要具有面性、均匀、随机分布的特点。

（3）整理原始数据统计表

收集因变量观测值及相应坐标，填入趋势分析原始数据统计表，其格式可参考表4.1 。

表4.1 趋势分析原始数据统计表格式

编 号

顺序号

原编号

控制点坐标

横坐标(u)

纵坐标(v)

变量

z

其中“控制点坐标’’可以是地理坐标，也可以是虚拟坐标系的相对坐标。 （4）建立原始数据矩阵

假设有n个控制点的地理坐标及地质变量观测值，则原始数据矩阵为：

?u1?u2???u?nv1v2?vnz1?z2? ??zn??（5）统一量纲

为了不改变地质变量的值，仅对控制点的坐标值进行以下变换（先中心化，在均匀化）：

ui??c?vi??c?1nnui?uuvi?vv?c?(?c?(uiuviv?1)(i?1,2,?,n) ?1)其中，u?i?1?ui，v?1ni?1?vi；c为常系数，其目的在于使ui?，vi?与zi具有相同的数量级。

n例如，在对间隔型（标高、水位等）进行趋势分析时，因为zi一般有三位整数，故将c取作100即

15

可；对比例型数据（煤厚、灰分、显微组分含量等）进行趋势分析时，可视zi的数量级，为c值在1~100中间取值。

若自选坐标系，控制点的坐标为相对坐标，可通过适当选取坐标系的单位，达到与地质变量统一量纲的目的，因而不必再施加上述变换。

2、趋势方程的参数估计

纪正规方程组为：

?1?u?i?X???vi????vni??ui2?ui?uivi??uivin?vi?uivi?v?2i?????vin?1n?vi??b0???zin??uivi?b1???uizi????n?1v?i?b2??vizi??????????nn?n?bvzi??l??i?v?i???? （4-1） ????解方程组可求出系数b0，b1，b2，?，bl，从而得到趋势面方程：

??i?b0?b1ui?b2vi???blvi （4-2）

n其中，l?0.5?(n?1)(n?2)

3、趋势与剩余、异常与“噪声”的分离

（1）由趋势方程计算各控制点的趋势值，再由观测值、趋势值求得各点的剩余值：

Ri?zi???i (i?1,2,?,n) （4-3）

（2）计算正剩余的平均值

R??1m??zj （4-4）

j?1M?或正剩余的二倍标准差

2S??21m?1?(Rj?R) （4-5）

j?1m??2（3）计算异常点和异常值

?选R?或2S?作为异常限，记为E。若正异常剩余与异常限之差ri?Rj?E?0，则第i点为正

异常点，ri为正异常值；若负剩余与异常限之和ri?Rj?E?0，则第i点为负异常点，ri为负异常值。

－4、趋势方程的显著性检验

（1）计算趋势平方和S趋、剩余平方和S剩、拟合度C及F检验统计量

S趋??(??i?z) （4-6）

i?1n2 16

S剩??(zi???i) （4-7）

i?1n2S总?S趋?S剩 （4-8）

C?F?S趋S总?100% （4-9）

S趋pS剩n?p?1 （4-10）

（2）给定置信水平?，查F临界值表的F?(p,n?p?1)。若F?F?(p,n?p?1)，则认为趋势方程在置信水平下是显著的。

5、绘制趋势图和偏差图，并进行地质解释 三、程序设计框图

要求：针对算例编写，不要求通用性，这样可以少占用一些时间。

建议以实验一的原始数据格式输入，首先要对数据进行预处理，处理完成后，形成实验一的原始数据格式的新文件（详见图1－1）。以下步骤如流程图（图4－1）。

开始打开标准格式的原始数据文件Nrj?R?j?E?0？Rj<0?Yrj?R－j?E?0？数据预处理和统计模块rj为正异常值输出变换后的标准格式数据文件rj为负异常值输出正负异常值?i?b0?b1ui?b2vi???blvin(i?0,1,?，n)求出z和?n?S?zi2?nz2?总?i?1n?22?S趋????i?nzi?1?计算?S剩?S总?S趋?C?SS（?100％）总趋?S趋p?F???S剩n?p?1?读入n对数据（ui，vi，yi）（i＝1，2，?，n）解正规方程组（5－1式）求得b0,b1,?,b（详细步骤见多元回归分析）l得方程??i?b0?b1ui?b2vi???blvin?i (i?1,2,?,n)计算剩余值：Ri?zi??计算正剩余的平均值：计算正剩余的二倍标准差：1M1m??2R????z?2S??2?(Rj?R)jmj?1m?1j?1输入F?（或者建立数据库，输 入?，系统自动查询F? ）F>Fα用户选择E＝R?或2S?Y输出C,F,S总,S趋,S剩,F?N结束

图4-1 趋势分析分析模块流程图

17

四、算例

（一）、在15*15km范围内，选择15个观测点测量某沉积物的粒径（坐标单位：km；粒径单位：mm）

坐标 粒径（z） 编号 粒径（z） 编号 横坐标(u) 纵坐标(v) 横坐标(u) 纵坐标(v) mm mm km km km km 1 2 3 1.9 9 8 11 1.3 2 2 10 2.3 10 10 8 1.2 11 3 2 13 1.1 11 13 1.4 4 4 1 2.6 12 12 3 1.7 5 5 8 2.2 13 12 6 1.8 6 5 14 1.8 14 12 10 1.2 7 7 3 3.5 15 15 13 1.0 8 7 6 3.1 16 18 16 1.4 要求进行趋势面分析，具体如下： （1）建立合适的面趋势面方程； ?和e； （2）求出各点的Zii2

坐标 （3）计算F值并进行检验； （4）绘趋势图和偏差图； （5）进行地质解释。

（二）、对下表中煤层底板标高数据（Z8,Z9）进行趋势面分析。要求： （1）建立合适的面趋势面方程； ?和e； （2）求出各点的Zii（3）计算F值并进行检验；

（4）绘一次趋势图和一次偏差图，比例尺为1：5000； （5）进行地质解释。 编号 1. 2. 3. 4. 5. 6. 钻孔号`11a` `12a` `13` `21` `22` `23` Xi Yi Z8 Z9 Z10 Z11 247.5 217.5 155.0 540.0 486.5 431.5 987.5 886.0 611.0 1125.0 902.5 673.0 -38.0 -60.0 -105.0 -13.5 -54.5 -95.5 -88.0 -110.0 -155.0 -63.5 -104.5 -145.5 -38.0 -60.00 -105.0 -13.5 -54.5 -95.50 -38.0 -60.0 -105.0 -13.5 -54.5 -95.5 18

7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. `24` `25` `31`， `32`， `33`， `34`， `35`， `41`， `42`， `43`， `44`， `45`， `51`， `52`， `53`， `54`， `61`， `62`， `63`， `64`， `65`， 367.5 282.5 777.5 717.0 656.5 591.0 559.0 1055.0 986.0 922.0 858.0 804.5 1276.0 1214.5 1154.5 1060.0 1575.0 1502.5 1456.5 1385.0 1329.5 424.0 87.5 1049.5 80.0 536.5 284.0 67.5 1121.5 849.0 576.0 341.0 119.0 975.0 722.0 490.0 128.0 1107.5 826.0 624.5 388.0 171.0 -141.5 -204.0 -26.5 -71.0 -119.0 -166.0 -207.0 -63.5 -113.5 -160.0 -203.5 -254.0 -90.0 -134.0 -176.0 -240.5 -67.0 -116.0 -151.0 -196.0 -234.9 -191.5 -245.0 -76.5 -121.0 -169.0 -216.0 -257.0 -13.5 -63.50 -110.0 -153.5 -195.0 -40.0 -84.0 -126.0 -190.5 -17.0 -66.00 -101.0 -146.0 -184.9 -141.5 -204.0 -76.5 -121.0 -169.0 -216.0 -257.0 -63.5 -113.5 -160.0 -203.5 -245.0 -140.0 -184.0 -226.0 -290.5 -117.0 -166.0 -201.0 -246.0 -284.9 -141.5 -204.0 -76.5 -121.0 -169.0 -216.0 -257.0 -63.5 -113.5 -160.0 -203.5 -245.0 -40.0 -84.0 -126.0 -190.5 -17.0 -66.00 -101.0 -146.0 -184.9 19

试验五 聚类分析

一、目的：

通过对聚类分析简单程序设计及算例计算，掌握聚类分析方法步骤及其数学原理。

二、方法概要

1、合理地选择变量和样品，构造原始数据矩阵

设有n个样品，m个变量，则原始数据矩阵?x??xij??n?m

其中，xij为第i个样品第j个变量的观测值（i=1,2,?,n;j=1,2,?,m）。

2、数据预测理（统一量纲变换）

通常，对服从正态分布的变量进行标准化变换，对服从对数正态分布的变量进行对数化变换，而变量不服从正态分布，或各变量服从不同的分布时，则采用正规化处理（详细处理过程见实验一）。

3、计算聚类统计量

一般在进行R型聚类时，计算相关系数；进行Q统计量型聚类时，计算相似系数或距离函数，其中以斜交距离函数的分类效果较好。事实上，各种方法各有所长，在计算中可选用多个统计量分别进行计算，然后结合实际情况，选择合适的方法。

（1）R型聚类采用相关系数或相似系数，按下式计算相关系数：

??xrjk?i?1nij?xj??xik?xk2n????xi?1nij?xj??xik?xkSjj??j,k?1,2,?,m? （5－1）

??xi?1nij?xj???i?1?xik?xk?2?Skk其中，i?1,2,?,n样品，xj,xk分别为第j个变量和第k个变量的样本均值，Sjj,Skk分别为第j个变量和第k个变量的样本方差。

当原始数据经过标准化数据后，由于xj?xk?0，Sjj?Skk?1，（6－1）式可简化为

rjk?1nij?xni?1?xik?i?1,2,?,n样品，j,k?1,2,?,m变量? （5－2）

（2）相似系数（S）

m?xSik?j?1mij?xkjm?i,k2?1,2,?,n样品，j?1,2,?,m变量? （5－3）

?xj?12xj??xkjj?1 20

（3）欧氏距离系数（D）

dik???xj?1mij?xkj?2m?i,k?1,2,?,n样品，j?1,2,?,m变量? （5－4）

（4）斜交距离系数（PD）

???xPDIK?j?1l?1mmij?xkj??xil?xkl???m2jl ?i,k?1,2,?,n样品，j?1,2,?,m变量?（5－5）

?jl：为第j个变量与第l个变量的相关系数。

（5）误差平方和增量（ΔE）

若P点群与q点群合并为t点群，误差平方和增量为ΔEpq,则

?Epq?nq?ngnp?nq??xj?1mjp?xjq?2初次12??xj?1mjp?xjq?2 （5－6）

4、划分类别（点群归并）

可用一次计算或逐步计算聚类法对变量或样品进行分类。两种方法比较而言，后者较前者较为合理一些，但计算量比前者约大一倍。

（1）对相关系数（R）和相似系数（S）,用逐步计算法或加权平均迭次推算法 ①修改变量数据

?为 Q型： i ，k合并，则xik??xiknixij?nkxkjni?nk?j?1,2,?,m? （5－7）

式中， ni,nk分别为第i,k点群中样品的个数R型：j,k合并，则

x?jk?njxij?nkxiknj?nk。

?i?1,2,?,n? （5－8）

式中，ni,nk分别为第j,k点群中变量的个数。 ②加权平均迭次推算法

若P点群与q点群合并为t点群，则t点群与其他r点的群的相似性系数为：

Str?npSpr?nqSqrnp?nq?r?1,2,?,m?1? （5－9）

（2）对于D、PD和 ΔE采用刷新方程迭次推算法

22①D、PD的刷新方程（含：Str代表D或PD）

21

若P,q合并为t，则t与其他r点群的St为：

Str?npnp?nq?Spr?nqnp?nq?Sqr?np?ng?np?nq?2Spq （5－10）

式中：np,nq分别为p,q点群中的样品个数。 ②误差平方和增量（ΔE ）的刷新方程。

若p,q点群合并为t点群，则t与r点群合并的ΔE为

?Etr?1nt?nr??np?nr??Epr??nq?nr???Eqr?nr??Epq （5－11） ?式中，nt,nr,np,nq分别为t, r, p, q点群中样品的个数。

5、整理联结表，绘制谱系图 6、对聚类结果进行解释 三、程序设计框图

要求：对算例设计程序，不要求程序具有通用性。 聚类分析程序框图（见图5－1）

开始打开标准格式的原始数据文件形成数组x[n,m]对数化变换对数化变换对数化变换求相关系数式（6－2）形成相关系数矩阵[R]求相似系数式（6－3）形成相似系数矩阵[Q]求欧氏距离系数式（6－4）形成欧氏距离系数矩阵[D]求斜交距离系数式（6－5）形成斜交距离系数矩阵[Pd] 一次聚类法：从相似性矩阵中选出最大的一对变量或样品，形成第一级分类，选出次大的一对变量或样品，形成第二级分类， ? ? ，到所有变量或样品聚合成一类为止。 逐步聚类法：从相似性矩阵中选出最大的一对变量或样品，例如rst?max{rjk}(j,k?1,2,?,m)(s?t),则将它们合成一类，并从nsxsj?ntxtj[x]矩阵中划去第t行，将xsj?(j?1,2,?,m)取代[x]矩阵的第s行ns?nt分类方案(1)分类方案(2)分类方案(h)分类方案(q)经与实际情况对比，假设分类方案(h)为最佳根据分类方案（h）整理联结表，绘制谱系图进行地质解释 图5-1 聚类分析模块流程图

22

四、算例

（一）、在广西合山朔河煤矿，采得四个煤样，每个煤层测得Ag,Scq和灰分中SiO2含量。

煤样 1 2 3 4 Ag?00? Scg?00? SiO2?00? 67.90 62.58 51.00 63.64 36.38 37.16 39.97 21.72 7.78 8.01 6.00 9.95 用正规化数据的相似系数Q型点群分析。 利用标准化数据的相关系数做R型点群分析。

（二）中南工业大学地质系在某次二年级“地质填图”实习中，师生测得4条断层的几何要素如下表。试研究这4条断层之间的异同程度和内在联系。

断层代号 A B C D 倾向（°） 100 110 10 20 倾角（°） 50 60 15 10 破碎带宽（m） 10 8 5 10 水平断距（m） 200 168 56 40 23

试验六 两类判别分析

一、目的：

通过对两类判别分析程序设计及完成算例，掌握两类判别分析基本原理和方法。

二、方法概要

1、整理原始数据，构成数据矩阵，并进行必要的预处理

设有n个样品，每个样品测得p个变量的数值，且在这n个样品中，已知有nA个样品属于A类，nB个样品属于B类，nV个样品的归属待定，但必然属于A类或B类。此时，原始数据矩阵如下：

?xA11?xA21????xAnA1?xB11?x[X]??B21???xBnA1?xv11?xv21?????xvnA1xA12xA22?xAnA2xB12xB22?xBnA2xv12xv22?xvnA2xA1p??xA2p??????xAnAp??xB1p??xB2p??????xBnAp??xv1p??xv2p??????xvnAp???????A类样品??????B类样品 （6－1） ??????待判样品???为消除量纲的影响，对[X]矩阵进行统一量纲变换，变换后的矩阵仍以[X]表示。

2、利用A类和B类样品的变量观测数据（经过统一量纲变换），建立判别函数

（1）计算组内平均值和组内方差

xAj?1nAnA?xk?1Akj；xBj?nA1nBnB?xk?12Bkj(j?1,2,?,r) （6－1）

2BKjSji?1nA?nB?2nB[?(xAKj?xAj)?k?1?(xk?1?xBj)] （6－3）

（2）计算各变量的I值

Ij?(xAj?xBj)/Sji (j=1,2,?,r) （6－4）

2（3）挑选变量

根据Ij的大小，从r个变量中，挑选出P个变量（如Ij?0.01）参加建立判别函数。

3、建立判别方程

（1）计算P个变量的组内方差和协方差

24

Sji?1nA?nB?2nAnB[?(xAKj?xAj)(xAKj?xAj)?k?1?(xk?1BKi?xBi)(xBKj?xBj)] (i,j=1,2,?,P) （6－5）

（2）解正规方程组

?S11C1?S12C2?????S1pCp?xA1?xB1??S21C1?S22C2?????S2pCp?xA2?xB2??......??Sp1C1?Sp2C2?????SppCp?xAp?xBp （6－6）

P求出C1，C2?..Cp,获得判别函数y??Cj?1jxj

4、计算分界值（y0），得出判别规则

（1）计算判别值、类平均值和总平均值

PPjyAB??Cj?1nAxAkj;yBk??Cj?1nBjxBkj?(k?1,2,...,nA) （k=1,2,?,nB）（6－7）

yA?1nA?k?1yAk;yB?1nB?k?1yBk;y?nAyA?nByBnA?nB （6－8）

（2）设：yA?yB

y0?yA?(yB?yA)?A?A??B （6－9）

（3）总结出判别规则

若有yA?y0?yB，则判别规则如下：

若待判样品的判别值y?y0，则将其归入A类； 若待判样品的判别值y?y0，则将其归入B类；

5、待判样品判别归类

P计算未知样品得判别值yx??Cj?1jxj，其中(j=1,2,?,P)

若待判样品的判别值y?y0，则将其归入A类； 若待判样品的判别值y?y0，则将其归入B类；

6、显著性检验

设：H0:yA?yB （1）计算F值

25

F?nA(yA?y)?nB(yB?y)?(nA?nB?2)nA22?(yk?1Ak?yA)?2nB （6－10）

?(yk?1Bk?yB)2在给定?下，查F分布表得F?(1,nA?nB?2)，如果F?F?(1,nA?nB?2)，则判别函数有效，否则，判别函数无效。 （2）计算回代正确判别率（r）

yAk?y0者为正确，共有m1个，rA?m1nA?100%

yBk?y0者为正确，共有m2个,rB?m2nB?100%

正确判别率：

r?m1?m2nA?nB?100% （6－11）

三、程序设计框图

开始读入nv对代判样品数据（xvi1,xvi2,?,xvip）(i?1,2,?,nv)打开标准格式的原始数据文件求出yj??Cx(j?1,2,?,n)jjvj?1p数据预处理和统计模块输出变换后的标准格式数据文件yA?yB判别规则：y?y0，归入A类；否则，归入B类判别规则：y?y0，归入B类；否则，归入A类读入nA对数据（xAi1,xAi2,?,xAip）(i?1,2,?,nA)读入nB对数据（xBi1,xBi2,?,xBip）(i?1,2,?,nB)计算组内平均值和组内方差（6－1式）计算F值（6－10式）输入F?（或者建立数据库，输 入?，系统自动查询F? ）计算各变量的I值（6－2式）用户设定输入挑选变量的界限值，挑选大于界限值的变量F>FαY计算P个变量的组内方差和协方差（6－5式根据方程y?解正规方程组（6－6式）求得C0,C1,?,C（详细步骤见多元回归分析）Pp?Cx计算njjj?1pA?nB个样品的yN根据判别规则判别样品的归类，并统计归类正确的数mA和mBmA?mB?100%nA?nB得方程y??Cxjj?1j计算正确率：r?计算分界值y（8，9式）06－7，输出以上计算的所有结果结束图6-1 二类判别模块流程图

26

四、算例

今获得有关内陆泥炭和滨海泥炭得锶（Sr）和钡（Ba）的化验数据如下： 样品号 锶（Sr） 1 0.0012 2 0.0030 滨海 3 0.0003 泥炭 4 0.0052 5 0.0002 1 0.0009 2 0.0007 内陆 泥炭 3 0.0025 4 0.028 1 0.001 未知 2 0.005 个体 3 0.0007 试用两类判别分析判断（未知个体）煤层的成因类型。 已知类 钡（Ba） 0.0001 0.0005 0.0002 0.0018 0.0002 0.0022 0.0013 0.0110 0.0060 0.0005 0.00065 0.0002 Sr / Ba 12 6 1.5 2.88 1 0.41 0.54 0.23 0.47 2.00 7.69 3.50 27

试验七 贝叶斯多类判别分析

一、目的：

通过对多类线性判别分析程序设计和完成算例，掌握贝叶斯多类判别模型原理和工作方法步骤。

二、方法概要 1、原始数据获取

设有G类母体，从每个母体中取得ng个样品，每个样品测得p个变量，则原始数据为：

xgkj（g ? 1,2,?,G ；k ? 1,2,?,ng；j ? 1,2,?,p ）

GN??ng?1g（总样品个数）

2、准备工作

（1）计算诸变量的类平均值和总平均值

xgj?1ngng?xk?1gkj xj?1GG?g?1Xgj （7－1）

（其中g ? 1,2,?,G ；k ? 1,2,?,ng；j ? 1,2,?,p） xj?1NGnggkj??xg?1k?2或xj??G1GXgj （7－2）

g?1（2）计算组内离差矩阵（W）和总离差矩阵（T）

W?(?ij)pxp;T?(tij)pxp （7－3）

Gnggki?ij???(xg?1k?1Gng?xgi)(xgkj?xgj) （7－4）

tij???(xg?1k?1gki?xi)(xgkj?xj) (i,j ? 1,2,?,p) （7－5）

（3）求 W，T矩阵的逆矩阵（W-1,T-1） 及行列式值

3、判别分类

（1）计算判别系数Cjg和Cog及先验概率qg

ppijCjg??(N?G)??i?1?xgj?(N?G)???ij?xgj （7－6）

i?1(g ? 1,2,?,G ；j ? 1,2,?,p)

28

Cog??12p?Ci?1jgxgj (g ? 1,2,?,G ) （7－7）

qg?ngN, (g ? 1,2,?,G ) （7－8）

（2）检验P个变量的判别效果 用威尔克斯准则U?体的能力。

1WT来检验H0:a1?a2???ag，即检验p个变量对于区分G个母

计算F?1?U1a??p(G?1) （7－9）

Ua其中：

?p2(G?1)2?4当p2?(G?1)2?5?0时 ?a??p2?(G?1)2?5 当p2?(G?1)2?5?0时

??1??ka?P(G?1)2P?G2?1

k?N?1?

在给定?下，查F分布表得F?(p(G?1),?)，如果F?F?(p(G?1),?)，则判别函数有效，否则，判别函数无效。

（3）计算未知个体X ? (x1,x2,?,xp)的判别值：

pyg(X)?lnpg?Cog??Cj?1jgxj （g ? 1,2,?,G ） （7－10）

（4）对未知个体类别分类

若yg(X)?max{yg(X)} (g ? 1,2,?,G )，则将X样品划归第g个母体Ag。

*4、计算后验概率

P{Ag/X}?eGy'g(X) （7－11）

y'g(X)?ek?1式中：y'g(X)?yg(X)?yg(X)

* 29

5、正确判别率估计

设对已知类型n个样品判别归类后，有m个样品归类正确，则正确判别率

r?mn?100% （7－12）

三、程序设计框图

开始打开标准格式的原始数据文件数据预处理和统计模块输出变换后的标准格式数据文件读入n对数据（xgk1,xgk2,?,xgkp）(g?1,2,?,G;k?1,2,?,ng;j?1,2,?,p)计算诸变量的类平均和总平均值（7－1，2式）计算组内离差矩阵（W）和总离差矩阵（T）（7－3，4，5式）计算判别系数Cjg和Cog以及先验概率qg(g?1,2,?,G;j?1,2,?,p)(7?6,7，8式）计算F值（7－9式）输入F?（或者建立数据库，输 入?，系统自动查询F? ）F>FαY计算未知样品X(x1,x2,?,xp)的判别值yg(X（)7－20式）(g?1,2,?,G)求出上述yg(X)(g?1,2,?,G)的最大值yg(X)，则归入g类计算后验概率（7－11式）?N根据7－10式计算n个已知类型样品的yg(X)根据最大值确定其类别，并统计归类正确的样品数mm?100%n计算正确率：r?输出以上计算的所有结果结束 图6-1 贝叶斯多类判别模块流程图

30

四、算例

（一）、某煤矿矿井开采A、B、C三个煤层，由于断层破坏，掘进巷道遇到一层没不知属于哪一层，从而影响掘进工作正常进行。试用判别分析解决该煤层对比问题。为此，取若干煤样，经过化验获得如下数据：

煤层 A 样品 1 2 3 4 1 2 3 4 1 2 3 4 1 S 7.83 7.58 8.51 8.31 4.73 5.12 5.78 6.17 6.56 7.78 7.28 7.32 4.56 Al2O5 23.35 23.20 23.89 24.00 38.92 37.84 37.19 37.40 14.59 15.84 14.95 13.94 39.02 C2O 2.74 3.15 4.19 4.32 2.40 3.23 3.36 3.30 2.33 4.39 9.32 3.33 1.36 已知类型 B C 未知个体 X 试建立多类线性判别函数，并对未知个体进行判别，求出X属于各母体的后验概率。

31

试验八 有序地质量最优分割法

一、目的：

通过对有序地质量最优分割法简单程序设计及完成算例，掌握该法的数字原理与方法。

二、方法概要 1、取得原始数据

设有N各有序样品，每个样品测得P各变量，则有原始数据矩阵

?XX?1112...X1P?X??X21X22...X?1P?????? ??XN1XN2...X?NP?2、数据正规化

ZXxj?min{Xij}?i?1,2,.N..?,ij?max{X ?ij}?min{Xij}?j?1,2,.P..? ?,3、计算段直径矩阵D

jp2d(i,j)?????Z???Z?(1,j)???1??? 1?i?j?N 1其中：

Z1j?(i,j)?j?i?1?Z??

??1得

?d(1,1)d(1,2)...d(1,N)??D??d(2,2)...d(2,N)??????NXN ??d(N,N)??

4、由D矩阵计算全部分段的组内离差平方和，并求出各段的最优分割。（1）最优二段分割

1） 由D矩阵对每个m = N,N-1，……，2计算相应的组内离差平方和。

Wm(2;j)?d(1,j)?d(j?1,m) (j=1,2,…,m-1)

32

8－1）8－2）8－3）8－4）8－5） （ （

（ （

（

2） 找出最优二段分割点

Wm(2;?1(m))?min{Wm(2;j)} （1?j?m?1） （8－6）

3） N个样品的最优二段分割为：

{x1,x2,…..xa1(N)} { xa1(N)+1…..xN} （8－7）

最优二段分割点为a1(N)。

（2）最优三段分割

1） 由D和最优二段分割结果，对每个m = N,N-1,…..3 计算相应的组内离差平方和。

Wm(3;?1(j),j)?Wj(2;?1(j))?d(j?1,m) （j=2,3,…,m-1） （8－8）

2） 确定最优三段分割及分割点

Wm(3;?1(m),?2(m))?min{Wm(3;?1(j),j)} （2?j?m?1） （8－9）

N个点的最优三段分割为：

{x1,x2,x?1(N)?1}{x?1(N)?1...x?2(N)}{x?(N)?1...xN} （8－10）

最优三段分割点为：?1(N),?2(N)

（3）最优K段分割

1） 由D和K－1段分割计算的结果，对每个m = N,N-1,….k分别计算相应的组内离差平方和。

Wm(k;?1(j),?2(j),...,?k?2(j),j)?Wj(k?1;?1(j),?2(j),...,?k?2(j))?d(j?1,N)

(j=k-1,k,…,m-1） （8－11）

2） 确定最优K段分割及分割点

Wm(k;?1(m),?2(m),...,?k?1(m))?min{Wm(k;?1(j),?2(j),...,?k?1(j),j)}

（k?1?j?m?1） （8－12）

N个样品的最优K段分割为：

{x1...x?1(N)}{x?1(N)?1...x?2(N)}...{x?(K?2)(N)...x?(K?1)(N)}{x?(K?1)...xN} （8－13）

最优K段分割点为：?1(N),?2(N),...?(k?1)(N)

5、给定K值或绘制W-K曲线，确定最终分段段数。

33

三、程序设计框图

开始挑选最小变差（8－5式）打开标准格式的原始数据文件（包括N和最高段数K）确定最优L段分割（8－6式）数据正规化（实验一）L＝L+1计算L段分割的组内离差平方和（8－11式）Y计算段径矩阵D数据正规化（实验一）挑选最小变差（8－12式）计算段径矩阵D（8－4式）确定最优L段分割（8－13式）L＝2计算L段分割的组内离差平方和（8－5式）L≤KN输出2至K段的最优分割点和最优分割结束

四、算例

(一)、某煤矿所采煤层牌号为主焦煤。在煤巷中见一火成岩墙侵入煤层，致使煤质发生变化，为弄清煤质变化情况，从火成岩附近，每隔0.5m依次取一煤样，现有6个有序煤样的镜煤最大反射率测定数据为

R0(%)??x1,x2,x3,x4,x5,x6???3.20,2.95,2.35,1.80,1.50,1.45?

试将这些有序数据进行最优三段分割，并说明岩墙对煤质有何影响？

34

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