Transcript 回归分析方法 - 数学建模学习资源网站

数据处理专题
一. 一元数据处理方法
二. 多元数据处理方法
三.如何写好建模竞赛论文
数据处理是指用简明而严格的方
法把获得的实验数据所代表的事
物内在的规律提炼出来，得出结
果的加工过程，包括数据记录、
描绘曲线，从带有误差的数据中
提取参数，验证和寻找经验规律
，外推实验数据等等。本章介绍
一些最基本的数据处理方法。
1. 一元数据处理方法
1.插值
2.拟合及线性回归
在解决实际问题的生产（或工程）实践和科
学实验过程中，通常需要通过研究某些变量之
间的函数关系来帮助我们认识事物的内在规律
和本质属性，而这些变量之间的未知函数关系
又常常隐含在从试验、观测得到的一组数据之
中。因此，能否根据一组试验观测数据找到变
量之间相对准确的函数关系就成为解决实际问
题的关键
例如在工程实践和科学实验中，常常需要从一组试验观测
数据(xi ,yi ) ，i = 0,1,....,n之中找到自变量x与因变量y 之间的
函数关系，一般可用一个近似函数y = f (x)来表示。函数
y = f (x)的产生办法因观测数据和要求不同而异，通常可采用
数据拟合与函数插值两种办法来实现。
数据拟合主要是考虑到观测数据受随机观测误差的影响，
进而寻求整体误差最小、能较好反映观测数据的近似函数y =
f (x)，此时并不要求所得到的近似函数y = f (x)满足
yi= f (xi) ， i = 0,1,…,n。
函数插值则要求近似函数y = f (x)在每一个观测点 xi 处一定
要满足y i= f (xi) ， i = 0,1,…,n ，在这种情况下，通常要求
观测数据相对比较准确，即不考虑观测误差的影响。
在实际问题中，通过观测数据能否正确揭示某些变量之
间的关系，进而正确认识事物的内在规律与本质属性，
往往取决于两方面因素。其一是观测数据的准确性或准
确程度，这是因为在获取观测数据的过程中一般存在随
机测量误差，导致所讨论的变量成为随机变量。其二是
对观测数据处理方法的选择，即到底是采用插值方法还
是用拟合方法，插值方法之中、拟合方法之中又选用哪
一种插值或拟合技巧来处理观测数据。插值问题忽略了
观测误差的影响，而拟合问题则考虑了观测误差的影响
。但由于观测数据客观上总是存在观测误差，而拟合函
数大多数情况下是通过经验公式获得的，因此要正确揭
示事物的内在规律，往往需要对大量的观测数据进行分
析，尤为重要的是进行统计分析。统计分析的方法有许
多，如方差分析、回归分析等。
数据拟合虽然较有效地克服了随机观测误差的影响，
但从数理统计的角度看，根据一个样本计算出来的拟合
函数（系数），只是拟合问题的一个点估计，还不能完
全说明其整体性质。因此，还应该对拟合函数作区间估
计或假设检验，如果置信区间太大或包含零点，则由计
算得到的拟合函数系数的估计值就毫无意义。这里所采
用的统计分析方法就是所谓的回归分析。另外还可用方
差分析的方法对模型的误差作定量分析。
对于插值方法，本文简单介绍最常用的插值法的基本
结论及其Matlab实现问题。由于数据拟合问题必须作区
间估计或假设检验，所以除了介绍最基本的数据拟合方
法——最小二乘法的基本结论及其Matlab实现问题外，
我们专门介绍了对数值拟合问题进行区间估计或假设检
验的统计方法。
即介绍回归分析方法及其Matlab实现。
数据处理问题通常情况下只是某个复杂实际问
题的一个方面或部分内容，因而这里所介绍的数
据处理方法——函数插值和数据拟合的方法（包
括回归分析）通常只能解决实际问题中的部分问
题——计算问题。一般来说，对实际问题进行数
学建模需要用到多方面知识，只有很少的情况下
可以单独使用本章所介绍的内容，故我们最后以
修改后的美国91年数学建模A题为例说明如何使用
数值计算知识建立数学模型，从而解决实际问题
的方法。
1、插 值 法
在生产和实验中，常常需要根据一张表格表示的函
数推算该表中没有的函数值.解决此类问题的简单途径之
一利用插值法。
插值在数学发展史上是一个老问题，它是和Gauss,
Lagrange, Newton等在著名数学家连在一起的。它最初
来源于天体计算——由若干观测值计算人一时刻星球的
位置。现在，插值法在工程技术和数据处理有许多直接
应用，而且也是数值积分、数值微分的基础。
1.1 插值概念与基础理论
1.1.1 插值问题的提法
对于给定的函数表
x
x0
Y=f(x) y0
x1
y1
……. xn
…….. yn
（1）
(其中 y  f ( x )在[a,b]上连续， x0， x1，…，xn 是 [a,b]上的
n+1个互异的点)，在某函数类{(x) }中求一个函数(x) ，使
(xi)=yi ,
(i=0,1,2,…,n)
并用函数(x) 作为函数 y=f(x) 的近似函数，即
y= f(x)  (x) , ( x∈[a,b] )
(2)
这类问题称为插值问题。 [a,b]称为插值区间， x0 , x1,
... , xn 称为插值节点，（2）称为插值条件，插值条件是
选择近似函数的标准，满足此条件的近似函数 (x) 称为
插值函数， f(x) 称为被插值函数。
函数类{(x) }有多种取法，常用的有代数多项式、
三角函数和有理函数。
最简单的插值函数是代数多项式，相应的插值
问题称为多项式插值。
求 f (x) 的n次插值多项式
通过曲线
y  f (x)
y  pn ( x )
的几何意义，就是
上的若干个节点，作一条代数曲线 y  pn ( x)
来近似代替曲线 y  f (x) 。如图所示。
§1.2 插值多项式的求法
在前面讨论插值多项式的存在唯一性时，实际上已提
供了它的一种求法，即通过求解线性方程组来确定其系数ai
(i=0,1,2,…,n)
但是这种方法不仅计算量大，而且因不能获得简明的
表达式而给理论和应用研究带来不便。在这里我们学习两
种简便而实用的求答。
1.2.1 拉格朗日插值多项式
在线性代数中知道，所有次数不超过n次的多项式构
成一个n+1维线性空间。其基有各种不同的取法。因此
尽管满足条件（4）的n次插值多项式是唯一的，然而它
的表达式可以有多种不同的形式。如果取满足条件：
0, i  k
（9）
l k( x i )  1, i  k
l x ,l1 x ,l 2 x ,,l n x  作为上述
的一组n次多项式 0
线性空间的基，则容易看出
l0 ( x ) y0  l1 ( x ) y1 
n
 l n ( x ) yn   yk l k ( x )
k 0
（10）
是一个次数不超过n的多项式。且满足插
值条件（4）。
因此，由n+1个代数多项式 l0 x,l1 x,l2 x,,ln x
线性生成的多项式（10）就是满足插值条件的n次插值多
项式。
满足条件（9）的多项式 l0 x,l1 x,l2 x,,ln x
称为n+1个节点的n次基本插值多项式（或n次基函数）
显然，求拉格朗日多项式的关键是求n次插值基函数。
因为
l k x  为n次多项式，且
0, i  k
l k( x i )  1, i  k
因此，可设
lk  x   Ak ( x  x0 )( x  x1 )...( x  xk 1 )( x  xk 1 )...( x  xn )
n
( x  x0 ) ( x  xk 1 )( x  xk 1 ) ( x  xn )
lk ( x ) 
( xk  x0 ) ( xk  xk 1 )( xk  xk 1 ) ( xk  xn )
两种特殊的Lagrange插值多项式
1.线性插值(两点插值)
最简单的插值是线性插值(此时n=1), 这时插值问题就是
求一次多项式
P1(x)=a0+a1x
使它满足条件
这时
l0 ( x ) 
x  x1
x0  x1
P1(x0)=y0 ,
l1 ( x ) 
P1(x1)=y1 ,
x  x0
x1  x0
于是线性插值多项式为
x  x0
x  x1
L1 ( x )  y0
 y1
x0  x1
x1  x0
y  y0
即
Ln ( x )  y0  1
( x  x0 )
x1  x0
它就是通过M0(x0,y0)和M1(x1,y1)两点的线段.
2.抛物插值
线性插值仅仅用两个节点以上的信息，精确度较差。为
了提高精确度，我们进一步考察以下三点的插值问题(n=2)：
这时 l ( x )  ( x  x1 )( x  x2 )
0
( x0  x1 )( x0  x2 )
( x  x0 )( x  x2 )
l1 ( x ) 
( x1  x0 )( x1  x2 )
( x  x0 )( x  x1 )
l2 ( x ) 
( x2  x0 )( x2  x1 )
由此得到抛物插值多项式 L2 ( x)  y0l0 ( x)  y1l1 ( x)  y2l2 ( x)
抛物插值又称三点插值.
例1 已知 y  ln x 的函数表
10
11
x
y
2.3026
2.3979
12
2.4849
13
2.5649
14
2.6391
分别用拉格朗日线性和抛物线插值求 ln(11.5)的近似值，
并估计误差。
%lagrange插值法的程序
function y=lagrange(x0,y0,x);
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
clear
x0=[10 11 12 13 14 ];
y0=[2.3026 2.3979,2.4849,2.5649 2.6391];
x=10:0.1:15;
y=lagrange(x0,y0,x);
plot(x0,y0,’+’,x,y)
1901年龙格(Runge) 给出一个例子:
f ( x) 
1
1  x2
定义在区间[-1，1]上，这是一个光滑函数，它的任意阶导
数都存在，对它在[-1，1]上作等距节点插值时，插值多项式情况，见图:
2
f(x)=1/(1+x2)
从图中，可
见，在靠近
-1或1时，
余项会随n
值增大而增
大，如
P12(0.96)=3
×6!但
f(0.96)=0.25
n=10
1.5
n=2
1
n=4
0.5
0
n=6
-0.5
n=8
-1
-1.5
-5
-4
-3
-2
-1
0
1
2
3
4
5
从图中，还可发现，在0附近插值效果是好的，即余
项较小，另一种现象是插值多项式随节点增多而振动更多。
这种插值多项式当节点增加时反而不能更好地接近被
插之数的现象，称为龙格现象。
上述现象和定理，告诉我们用高次插值多项式是不妥当
的，从数值计算上可解释为高次插值多项式的计算会带来舍
入误差的增大，从而引起计算失真。那么如何提高插值精度
呢？采用分段插值是一种办法。实践上作插值时一般只用一
次、二次最多用三次插值多项式。
1.2.2 分段线性插值
分段线性插值：
matalb调用格式：
yi=interp1(x,y,xi,’linear’)
分段线性插值的构造：
设f(x)是定义在[a,b]上的函数，在[a,b]上节点
a= x0< x1<x2<…<xn-1<xn=b,
的函数值为 y0 , y1 ,y2 ,…yn-1 ,yn 。
x,y为插值节
点，xi为待
求节点
(x)在每个子区间[xi , xi+1](i=0,1,2,,n-1)上是一次插
值多项式;
x  xi 1
x  xi
 ( x )  yi
 yi 1
xi  xi 1
xi 1  xi
xi ,  x  xi 1
这种分段低次插值称为分段线性插值.在几何上就是用
折线段带代替曲线,故分段线性插值又称为折线插值.
实际上是连接点( xk , yk ) , i  0,1,
, n的一条折线
分段线性插值曲线图：
1
0.8
0.6
曲线的光滑性较差
0.4
0.2
在节点处有尖点
0
-0.2
但如果增加节点的数量
-0.4
-0.6
减小步长,会改善插值效果
-0.8
-1
-4
-3
-2
-1
0
1
2
3
4
例1 已知 y  ln x 的函数表
10
11
x
y
2.3026
2.3979
12
2.4849
13
2.5649
14
2.6391
分别用拉格朗日线性和抛物线插值求 ln(11.5)的近似值，
并估计误差。
clear
x0=[10 11 12 13 14 ];
y0=[2.3026 2.3979,2.4849 2.5649 2.6391 ];
x=10:0.1:15;
y1=interp1(x0,y0,x,’linear’);
yy1=interp1(x0,y0,11.5,’linear’);
y2=interp1(x0,y0,x,‘cubic');
yy2=interp1(x0,y0,11.5,‘cubic');
subplot(1,2,1)
plot(x0,y0,'+',x,y1,11.5,yy1,’rO’)
title('Piecewise linear')
subplot(1,2,2)
plot(x0,y0,'+',x,y2,11.5,yy2,’rO’)
title('Piecewise cubic')
分段二次插值即：选取跟节点x最近的三个节点
xi-1,xi, xi+1进行二次插值，即在区间[xi-1, xi+1]，取：

i 1 ( x  x ) 
i
j 

f ( x )  L2 ( x )   yk 
 j  i 1 ( x k  x i ) 
k  i 1
 j  k

i 1
这种分段的低次插值叫分段二次插值，在几何上
就是用分段抛物线代替y=f(x)，故分段二次插值
又和分段抛物插值。
matlab调用格式
yi=interp1(x,y,xi,’cubic’) %二次多项式插值
1.3 三次样条插值
什么是样条: 是 指飞机或轮船等的制造过程中为描绘
出光滑的外形曲线(放样)所用的工具
样条本质上是一段一段的三次多项式拼合而成的曲线
在拼接处,不仅函数是连续的,且一阶和二阶导数也是连续的
1946年,Schoenberg将样条引入数学,即所谓的样条函数
1.4.1、三次样条插值函数
定义1.
a  x0 , x1 ,, xn  b为区间[a, b]的一个分割
如果函数S( x)在区间[a, b]上满足条件:
(1) S( x), S( x), S( x)都在区间[a, b]上连续,即S( x)  C 2 [a, b]
(2) S( x)在每个小区间[ xk , xk  1 ]上都是三次多项式
(3) 如果函数f ( x)在节点x0 , x1 ,, xn处的函数值为
f ( x j )  y j , j  0,1,, n
而S ( x)满足 S( x j )  y j , j  0,1,, n
------(1)
则称S( x)为f ( x)在[a, b]上的三次样条插值函数
由条件（2），不妨将s( x)记为
s( x) ={ si( x) , x [ xi , xi1 ], i  1,2,..., n}
si ( x)  ai x  bi x  ci x  di
3
2
(1)
其中ai，bi，ci，di为待定系数，共4n个。
由条件（1），有
si ( xi )  si +1 ( xi )

si ( xi )  si +1 ( xi )
s( x )  s ( x )
i +1
i
 i i
(i  1, 2,..., n  1)
(2)
由条件（3），有
si ( xi )  yi
(i  0,1, 2,..., n)
(3)
容易看出，（2）和（3）共有4n  2个方程，为
确定s( x) 的4n个待定参数，尚需再给出两个条件，
即所谓边界条件。通常使用的边界条件有以下三类：
第一类边界条件是
s( x0 ) =f 0

s( xn ) =f n
f 0, f n为给定的值。
第二类边界条件是
s( x0 ) =f 0

s( xn ) =f n
f 0，f n为给定的值。
当s( x0 ) =s( xn )=0时，样条函数在两端点不受力，呈
自然状态，故称之为自然边界条件。
第三类边界条件是周期条件。设f ( x) 是周期函数，
不妨设以xn - x0为一个周期，这时s ( x )也应以xn - x0
为周期的周期函数，于是 s( x) 在端点处满足条件：
s( x0 +0) =s( xn - 0) ; s( x0 +0) =s( xn - 0) ; s( x0 +0) =s( xn - 0) .
利用4n个条件求出三次样条函数s( x)的4n个待定常数，
直接求解计算量很大，通常利用Mat l ab软件求s( x).
例2　设f ( x) 为定义在[ 0, 3] 上的函数，有下列函数值表
xi
yi
0
0
1
0.5
2
2
3
1.5
且f ( x0 )  0.2, f ( x3 )  1, 试求区间[0,3]上满足上述
条件的三样次条插值函数s( x).
解 mat l ab求解s( x) . 程序为：
clear
x0=[0 1 2 3];y0=[0 0.5 2 1.5];x=0:0.1:3;
pp1=csape(x0,y0,’complete’);
y3=ppval(pp1,x);%计算插值函数在x处的值
plot(x0,y0,’+’,x,y3,’r’)
pp=csape( x, y, ' compl et e' )
[ br eaks, coef s, npol ys, ncoef s, di m] =unmkpp( pp)
上述程序中函数csape( x, y, ' complete ')是指定边界条件的
样条插值函数，csape返回一个包含三次样条插值的pp形，
或者说是分段多项式的结构。这个结构包含了计算用户希
望的任何插值点数值的三次样条值需要的所有信息。
字符串' compl et e' 表示所给边界条件是第一类边界条件；
若将' compl et e' 换成' second' 表示第二类边界条件；
' per i odi c' 表示第三类边界条件。
在第一类边界条件和第二类边界条件时，边界条件值放
在y的第一个分量和最后一个分量的位置上。
这样y的分量的个数比x的分量的个数多２。例如，在
例１中，
x=[ 0 1 2 3]
y=[ 0. 2 0 0. 5 2. 0 1. 5 - 1]
周期边界条件时，无需指定边界条件值。
在计算一个三次样条表达式的时候，必须将pp形中
不同域提取出来进行计算，这个过程可以由函数unmkpp
完成，该函数的使用方法为：
[ br eaks, coef s, npol ys, ncoef s, di m] =unmkpp( pp)
其中输入变量pp是样条插值函数csape的输出变量，
unmkpp的输出变量有５个：
br eaks, coef s, npol ys, ncoef s, di m。
br eaks：包含了插值节点；
coef s：是一个矩阵，其第i 行是第i 个多项式的系数：a0 , a1 , a2 , a3 .即
3
2
s　(
x
)

a
(
x

x
)

a
(
x

x
)
 a2 ( x  xi 1 )  a3
i
0
i 1
1
i 1
中的系数。
npolys : 是多项式的个数；
ncoefs : 是每个多项式系数的个数；
dim : 是样条的维数。
其中输入变量pp是样条插值函数csape的输出变量，
unmkpp的输出变量有５个：
br eaks, coef s, npol ys, ncoef s, di m。
br eaks：包含了插值节点；
coef s：是一个矩阵，其第i 行是第i 个多项式的系数：a0 , a1 , a2 , a3 .即
3
2
s　(
x
)

a
(
x

x
)

a
(
x

x
)
 a2 ( x  xi 1 )  a3
i
0
i 1
1
i 1
中的系数。
npolys : 是多项式的个数；
ncoefs : 是每个多项式系数的个数；
dim : 是样条的维数。
运行结果如下：
pp=
f or m: ' pp'
br eaks: [ 0 1 2 3]
coef s: [ 3  4 doubl e]
pi eces: 3
or der : 4
di m: 1
br eaks=
npol ys=
0 1 2 3
coef s=
0. 2800 - 0. 1800 0. 2000
- 1. 0400
3
0
1. 2600 1. 2800 0. 5000
0. 6800 - 1. 8600 0. 6800 2. 0000
ncoef s=
4
di m=
1
因此所求的三次样条函数为
0. 28x3  0.18 x 2  0.2 x,

s( x) = 1.04( x  1)3  1.26( x  1) 2  1.28( x  1)  0.5,
0. 68( x  2)3  1.86( x  2) 2  0.68( x  2)  2.0,

x  [0,1]
x [1, 2]
x  [2,3]
例２　已知函数值表
xi
yi
1
1
2
3
4
4
5
2
试求在区间１,
[ 5]上满足上述函数表所给出的插值条件的
三次样自然样条插值函数s( x).
解　用Mat l ab求解s( x) . 程序为：
x=[ １　２　４　５] ;
y=[ ０　１　３　４　２　０] ;
pp=csape( x, y, ' second' )
[ br eaks, coef s, npol ys, ncoef s, di m] =unmkpp( pp)
运行结果如下：
pp=
f or m: ' pp'
br eaks: [ 1 2 4 5]
coef s: [ 3  4 doubl e]
pi eces: 3
or der : 4
di m: 1
br eaks=
1 2 4 5
coef s=
- 0. 1250
0
2. 1250
1. 0000
- 0. 1250 - 0. 3750 1. 7500
3. 0000
0. 3750 - 1. 1250 - 1. 2500
4. 0000
因此所求的三次样条函数为
- 0. 125( x  1)3  2.125( x  1)  1,

s( x) = 0.125( x  2)3  0.375( x  2)2  1.75( x  2)  3,
0. 375( x  4)3  1.125( x  4)2  1.25( x  4)  4,

x  [1, 2]
x  [2, 4]
x  [4,5]
一维插值总结
插值函数一般是已知函数的线性组合或者称为加权平均。在已知数据点较少
时，插值技术在工程实践和科学实验中有着广泛而又十分重要的应用。例如在
信息技术中的图像重建、图像放大过程中为避免图像失真、扭曲而增加的插值
补点，建筑工程的外观设计，化学工程试验数据与模型分析，天文观测数据、
地理信息数据的处理，社会经济现象的统计分析等方面，插值技术的应用是不
可或缺的。
插值技术（或方法）远不止这里所介绍的这些，但在解决实际问题时，对于一
位插值问题而言，前面介绍的插值方法已经足够了。剩下的问题关键在于什么
情况下使用、怎样使用和使用何种插值方法的选择上。
拉格朗日插值函数在整个插值区间上有统一的解析表达式，其形式关于节点对
称，光滑性好。但缺点同样明显，这主要体现在高次插值收敛性差（龙格现
象）；增加节点时前期计算作废，导致计算量大；一个节点函数值的微小变化
（观测误差存在）将导致整个区间上插值函数都发生改变，因而稳定性差等几
个方面。因此拉格朗日插值法多用于理论分析，在采用拉格朗日插值方法进行
插值计算时通常选取n < 7。分段线性插值函数（仅连续）与三次样条插值函数
（二阶导数连续）虽然光滑性差，但他们都克服了拉格朗日插值函数的缺点，
不仅收敛性、稳定性强，而且方法简单实用，计算量小。因而应用十分广泛。
2、数 据 拟 合
在科学计算中经常要建立实验数据的数学模型。给定函数的实
验数据，需要用比较简单和合适的函数来逼近（或拟合）实验
数据。这种逼近的特点是：
(a) 适度的精度是需要的；
(b) 实验数据有小的误差；
(c) 对于某些问题，可能有某些特殊的信息能够用来选择实验数
据的数学模型。
逼近离散数据的基本方法就是曲线拟合，常采用最小二乘拟合
曲线拟合问题的数学描述是，已知一组（二维）数据(xi,yi ) ，
i = 1,2,。。。,n（即平面上的n个点(xi, yi ) ，i = 1,2,。。,n）， x i
互不相同。寻求一个函数（曲线） y = f (x)，使f (x)在某种准则下
与所有数据点最为接近，即曲线拟合得最好。
最小二乘拟合分为线性最小二乘拟合和非线性最小二乘拟合。
2.1线性最小二乘拟合（多项式拟合）方法
在线性最小二乘拟合中，用的较多的是多项式拟合。如果取
{ r1( x), ‥, rm+1( x)} ={1, ‥ ,xm } ，即用m 次多项式拟合给定
数据，则Matlab中有现成的函数
a=polyfit(x0,y0,m)，
其中输入参数x0,y0为要拟合的数据，m为拟合多项式的次数
，输出参数a为拟合多项式
y=amxm+…+a1x+a0系数a=[ am, …, a1, a0]。
多项式在x处的值y可用下面的函数计算
y=polyval(a,x)。
例4 某乡镇企业1990-1996年的生产利润如下表：
年份
1990 1991 1992 1993 1994 1995 1996
利润（万元） 70 122 144 152 174 196 202
试预测1997年和1998年的利润。
解 作已知数据的的散点图，
x0=[1990 1991 1992 1993 1994 1995 1996];
y0=[70 122 144 152 174 196 202];
plot(x0,y0,'*')
发现该乡镇企业的年生产利润几乎直线上升。因此，
我们可以用y = a1 x + a0 作为拟合函
数来预测该乡镇企业未来的年利润。编写程序如下：
x0=[1990 1991 1992 1993 1994 1995 1996];
y0=[70 122 144 152 174 196 202];
a=polyfit(x0,y0,1)
y97=polyval(a,1997)
y98=polyval(a,1998)
求得a1 = 20 ，a0 = -4.0705×104 ，1997年的生产利
润y97=233.4286，1998年的生产利润y98=253.9286
。
2 非线性最小二乘拟合
Matlab的优化工具箱中提供了两个求非线性最小二乘拟合的函数：
curvefit和leastsq。使用这两个命令时，都要先建立M文件fun.m，
但它们定义f (x)的方式是不同的。
1 curvefit
设已知xdata=(xdata1,xdata2,…,xdatan ),ydata=(ydata1,ydata2,…,ydatan ),
curvefit用以求含参量x（向量）的向量值函数
F(x,xdata)=(f(x,data1), …,f(x,xdata n )) T中的参变量x（向量），使得
Sum(F(x,xdatai)-ydatai)2最小
输入格式为：
（1）x=curvefit('fun',x0,xdata,ydata)；
（2）x=curvefit('fun',x0,xdata,ydata,options)；
（3）x=curvefit('fun',x0,xdata,ydata,options, 'grad')；
（4）[x，options]=curvefit('fun',x0,xdata,ydata,…)；
（5）[x，options,funval]=curvefit('fun',x0,xdata,ydata,…)；
（6）[x，options,funval,Jacob]=curvefit('fun',x0,xdata,ydata,…).
输出目标函数值格式：f=fun(x,xdata).
其中x0为迭代初值，options为控制参数。
3 非线性最小二乘拟合
Matlab的优化工具箱中提供了两个求非线性最小二乘拟合的函数：
curvefit和leastsq。使用这两个命令时，都要先建立M文件fun.m，
但它们定义f (x)的方式是不同的。
1 curvefit
设已知xdata=(xdata1,xdata2,…,xdatan ),ydata=(ydata1,ydata2,…,ydatan ),
curvefit用以求含参量x（向量）的向量值函数
F(x,xdata)=(f(x,data1), …,f(x,xdata n )) T中的参变量x（向量），使得
Sum(F(x,xdatai)-ydatai)2最小
输入格式为：
（1）x=curvefit('fun',x0,xdata,ydata)；
（2）x=curvefit('fun',x0,xdata,ydata,options)；
（3）x=curvefit('fun',x0,xdata,ydata,options, 'grad')；
输出目标函数值格式：f=fun(x,xdata).
其中x0为迭代初值，options为控制参数。
2 leastsq
设已xdata=(xdata1,xdata2,…,xdatan ),ydata=(ydata1,ydata2,…,ydatan ),
leastsq 用以求含参量x（向量）的向量值函数
输入格式为：
（1）x= leastsq ('fun',x0,options)；
（2）x= leastsq ('fun',x0,options, 'grad')；
例5 用下面一组数据拟合函数c(t) = a + be-0.02kt
中的参数a,b, k 。
t
cj×103
100 200 300 400 500 600 700 800 900 1000
4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59
1 用命令curvefit。此时
F(x，tdata)=(a+b e-0.02kt1，…，a+be-0.02kt10)T，x=(a，b，k)
（1） 编写M文件curvefun1.m
function f=curvefun1(x,tdata)
f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中x(1)=a；x(2)=b；x(3)=k；
（2） 输入命令
tdata=100:100:1000
cdata=1e03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];
x0=[0.2,0.05,0.005];
x=curvefit(‘curvefun1’,x0,tdata,cdata)
f=curvefun1(x,tdata)
即拟合得a=0.0070，b=-0.0030，k=0.0066
2.2.回归分析方法
2.1.1 引言
确定性关系
变量之间的关系
相 关 关 系
确定性关系
S  πr 2
相关关系
身高和体重
相关关系的特征是:变量之间的关系很难用一
种精确的方法表示出来.
确定性关系和相关关系的联系
由于存在测量误差等原因,确定性关系在实际
问题中往往通过相关关系表示出来;另一方面,当对
事物内部规律了解得更加深刻时,相关关系也有可
能转化为确定性关系.
回归分析——处理变量之间的相关关系的一
种数学方法,它是最常用的数理统计方法.
回归分析的任务——根据试验数据估计回归
函数;讨论回归函数中参数的点估计、区间估计;
对回归函数中的参数或者回归函数本身进行假设
检验;利用回归函数进行预测与控制等等.
回归分析方法
一元线性回归
数
学
模
型
及
定
义
模
型
参
数
估
计
检
验
、
预
测
与
控
制
性可
回线
归性
（化
曲的
线一
回元
归非
）线
多元线性回归
数
学
模
型
及
定
义
模
型
参
数
估
计
检多
验元
与线
预性
测回
归
中
的
逐
步
回
归
分
析
2.2一元线性回归分析
2.2.1 数学模型
例1 测16名成年女子的身高与腿长所得数据如下：
身高
143
145
146
147
149
150
153
154
155
156
157
158
159
160
162
164
腿长
88
85
88
91
92
93
93
95
96
98
97
96
98
99
100
102
以身高x为横坐标，以腿长y为纵坐标将这些数据点（xi，yi）
在平面直角坐标系上标出.
102
100
98
y   0  1 x  
96
94
92
90
88
散点图
86
84
140
145
150
155
160
165
一般地，称由 y   0  1 x   确定的模型为一元线性回归模型，
记为
 y   0  1 x  

2
E


0
,
D




（1）
固定的未知参数  0 、 1 称为回归系数，自变量 x 也称为回归变量.
Y   0  1 x ，称为 y 对 x 的回归直线方程.
一元线性回归分析的主要任务是：
1、用试验值（样本值）对  0 、  1 和  作点估计；
2、对回归系数  0 、  1 作假设检验；
3.对x  x0处对y做预测和区间估计
2.2.2 模型参数估计
（1）回归系数的最小二乘估计
有 n 组独立观测值，
（x1，y1）
，
（x2，y2）
，…，
（xn，yn）
 yi   0   x i   i , i  1, 2,..., n
设 
2
E


0
,
D



且 1 2，
...,  n 相互独立
 i
i
n
记
n
Q  Q(  0 ,  1 )       y i   0   1 x i 
i 1
2
i
i 1
最小二乘法就是选择  0 和 1 的估计 ˆ 0 ， ˆ1 使得
Q( ˆ 0 , ˆ1 )  min Q(  0 , 1 )
 0 , 1
2
为此，将上式分别对0，1求偏导，得
n
 Q
   2 ( yi   0  1 xi )
i 1
 0

n

Q

 2 xi ( yi   0  1 xi )
 1
i 1
令上两式为0，解得0，1，得
 ˆ 0  y  ˆ1 x

 ˆ xy  x y
 1  2
x  x2

n
ˆ1 
或
  x  x  y  y 
i 1
i
i
n
x  x 
i 1
2
i
其中
1 n
1 n
x   xi , y   y i
n i 1
n i 1
1 n 2
1 n
x   xi , xy   xi yi
n i 1
n i 1
2
用这种方法求出的估计ˆi (i  0,1)称为i的最小二乘
估计，简称LS 估计。（经验）回归方程为
yˆ  ˆ0  ˆ1 x  y  ˆ1 ( x  x)
显然，ˆ1是拟合直线的斜率，ˆ0是拟合直线的截距。n
个点( xi , yi )(i  1, 2,..., n)的几何重心( x , y )落在拟合直线上。
为了便于记忆，引入下列记号：
n
n

2
2
L

(
x

x
)
L

(
y

y
)

yy
i
 xx  i

i 1
i 1

n
 L  ( x  x )( y  y )

xy
i
i

i 1

（2）
这时ˆ1可简记为：
ˆ
1  Lxy Lxx .
注意：ˆ1
N ( 1 ,
2
n
2
(
x

x
)
 i
), 所以它是1的无偏估计，
i 1
同样，ˆ0也是 0的无偏估计。
（2）  2 的无偏估计
n
记

Qe  Q( ˆ 0 , ˆ1 )   yi  ˆ 0  ˆ1 xi
i 1
  (y
n
2
i 1
i
 yˆ i ) 2
称Qe为残差平方和或剩余平方和. 可以证明：
E(Qe )  (n  2)
于是 的无偏估计为
2
2
ˆ e2  Qe (n  2)
称ˆ e2为剩余方差（残差的方差），它是 2的无偏估计，
且ˆ 2分别与ˆ ，ˆ 独立。ˆ 也叫剩余标准差。显然，ˆ
e
0
1
e
越接近于0，说明线性回归方程（1）越显著。
e
2.2.3 回归方程的显著性检验
在实际工作中，实现我们并不能断定y与x之间
有线性关系，式( 1）只是一种假设。当然这个假设
不是没有根据，我们可以通过专业知识或散点图做
粗略判断。但在求出回归方程之后，还须对这种线
性回归方程同实际观测数据拟合的效果进行检验。
由（1）知，当| 1| 越大，y与x的变化趋势越明显；
反之，当当| 1| 越小，y与x的变化趋势越不明显；特
别当1 =0，y与x不存在线性关系。
对回归方程 Y   0  1 x 的显著性检验，归结为对假设
H0 : 1  0; 或 H1 : 1  0
进行检验.
假设 H 0 : 1  0 被拒绝，则回归显著，认为 y 与 x 存
在线性关系，所求的线性回归方程有意义；否则回归
不显著，y 与 x 的关系不能用一元线性回归模型来描
述，所得的回归方程也无意义.此时，可能：
（）
1 x对y没有显著影响，此时应去掉预报变量x;
( 2) x对y有显著影响，但这种影响不能用线性关系表示；
（3）除了x以外，还有其他不可忽略的变量对y的影响，
从而削弱了x对y的影响。此时应用多元回归模型。
设
n
n
Lyy   ( yi  y )   y i  ny
2
i 1
2
2
i 1
n
n
2
2
ˆ
ˆ
ˆ
Qe   ( yi  yi )   ( yi   0   1 xi )
i 1
n
剩余平方和
i 1
n
2
2
ˆ
ˆ
ˆ
U   ( yi  y )   (  0   1 x i  y )
i 1
总离差平方和
回归平方和
i 1
由于Lyy  Qe  U，因此U 越大，Qe 就越小。于是 U Lyy ，
越接近1，或 U Qe 越大，则线性回归效果越显著。
可以用三种不同方法进行检验，它们的本质
是相同的．这里介绍
F检验法
当 H 0 成立时，
U
~F（1，n-2）
F
Qe /(n  2)
n
2
ˆ
其中 U    yi  y  （回归平方和）
i 1
故 F> F1 (1, n  2) ，拒绝 H 0 ，否则就接受 H 0 .
回归系数的置信区间
 0 和 1 置信水平为 1-α的置信区间分别为
2
2 

1
x
1
x

, ˆ 0  t  (n  2)ˆ e

 ˆ 0  t  (n  2)ˆ e

1
1

n Lxx
n Lxx 

2
2
ˆ

ˆ
和   1  t  ( n  2)ˆ e / L xx ,  1  t  ( n  2)ˆ e / L xx 
1
1
2
2


 的置信区间为
 2 的置信水平为 1

Qe
Qe


  2 ( n  2) ,  2 ( n  2) 

 1 

2
 2

线性回归matlab
y  0  1x1
1．确定回归系数的点估计值：
b=regress( Y,
 ˆ0 
b   
ˆ

 1
 Y1 
Y 
Y   2
 
 
Yn 
X )
 1 x11 


x
1
21 
X 




 1 xn1 
b为一元函数y系数的估计值
2．求回归系数的点估计和区间估计、并检验回归模型：
[b, bint,r,rint,stats]=regress(Y,X,alpha)
回
归
系
数
的
区
间
估
计
残
差
显著性水平
（缺省时为0.05）
置
信
区
间
用于检验回归模型的统计量，
有三个数值：相关系数r 2、
F值、与F 对应的概率p
相关系数 r2 越接近 1，说明回归方程越显著；
F > F1-α（k，n-k-1）时拒绝 H0，F 越大，说明回归方程越显著；
与 F 对应的概率 p   时拒绝 H0，回归模型成立.
3．画出残差及其置信区间：
rcoplot（r，rint）
例1 测16名成年女子的身高与腿长所得数据如下：
身高
143
145
146
147
149
150
153
154
155
156
157
158
159
160
162
164
腿长
88
85
88
91
92
93
93
95
96
98
97
96
98
99
100
102
y   0  1 x  
clear
x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164]';
X=[ones(16,1) x];
Y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]';
[b,bint,r,rint,stats]=regress(Y,X);
b,bint,stats
rcoplot(r,rint)
figure(2)
a
z=b(1)+b(2)*x；
plot(x,Y,'k+',x,z,'r')
F=finv(0.95,1,14)%及求F0.05(1,16-2)
F (1, n  2)
b=
y  16.0730 0.7194x
-16.0730
0.7194
bint =
-33.7071 1.5612
0.6047 0.8340
- 33.7071 ˆ0  1.5612
0.6047  ˆ1  0.8340
相关系数 r2 越接近 1，说明回归方程越显著；
stats =
F > F1-α（k，n-k-1）时拒绝 H0，F 越大，说明回归方程越显著；
0.9282 180.9531
0.0000
与 F 对应的概率 p   时拒绝 H0，回归模型成立.
F0.05(1,16-2)=4.6
1.2.4 预测与控制
（1）预测:对固定的x值预测相应的y值
用 y0 的回归值 yˆ 0  ˆ 0  ˆ1 x0 作为 y0 的预测值.
y0 的置信水平为 1   的预测区间为
yˆ 0   ( x0 ), yˆ 0   ( x0 )
1 x  x 
其中  ( x0 )  ˆ e t  (n  2) 1   0
1
n
Lxx
2
2
特别，当 n 很大且 x0 在 x 附近取值时,
y 的置信水平为 1   的预测区间近似为


 yˆ  ˆ e u1 , yˆ  ˆ e u1 
2
2 

例1 测16名成年女子的身高与腿长所得数据如下：
身高
143
145
146
147
149
150
153
154
155
156
157
158
159
160
162
164
腿长
88
85
88
91
92
93
93
95
96
98
97
96
98
99
100
102
请你预测一下身高为170cm时该女子的腿长应该是在什么范围(置信度为0.05)
y   0  1 x  
y0 的置信水平为 1   的预测区间为
yˆ 0   ( x0 ), yˆ 0   ( x0 )
ˆ  Qe (n  2)
2
e
y  16.0730 0.7194x
1 x  x 
其中  ( x0 )  ˆ e t  (n  2) 1   0
1
n
Lxx
2
2
n
Lxx   ( xi  x )2
i 1
n

Qe  Q( ˆ 0 , ˆ1 )   yi  ˆ 0  ˆ1 xi
clear
x0=170
x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164];
y0=-16.073+0.7194*x0;
t=tinv(0.975,16-2);
Lxx=sum((x-mean(x)).^2)
d=sum(r.^2)/14*t*sqrt(1+1/16+(x0-mean(x))/Lxx)
yy=[y0-d,y0+d]
i 1
   ( y  yˆ )
2
n
i 1
i
i
2
clear
x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164]';
X=[ones(16,1) x];
Y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]';
[b,bint,r,rint,stats]=regress(Y,X);
x0=170
x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164];
y0=-16.073+0.7194*x0;
t=tinv(0.975,16-2);
Lxx=sum((x-mean(x)).^2)
d=sum(r.^2)/14*t*sqrt(1+1/16+(x0-mean(x))/Lxx)
yy=[y0-d,y0+d]
（2）控制
若要y   0  1 x  的值以1  的概率落在区间
（y, y) 之内，那么预报变量x应控制在什么范围内。
也就是说要求出区间( x, x) , 使当x ( x, x) 时，对
应的y值以1  的概率落在区间（y, y) 之内。这就
是控制问题。控制问题可以看成是预测问题的反问题。
只要控制 x 满足以下两个不等式
yˆ   ( x)  y , yˆ   ( x)  y 
要求 y   y   2 ( x) .若 yˆ   ( x)  y , yˆ   ( x)  y  分别有解 x 
和 x  ，即 yˆ   ( x)  y , yˆ   ( x)  y  .则 x , x  就是所求的 x
的控制区间.
一般说来，要解出x和x很复杂。若样本容量很
大，且x在x附近，问题可以大大简化。这时我们可以
利用预测区间[ yˆ  ˆ eu  , yˆ  ˆ eu  ]来进行控制。
1
2
1
2
1.3 可线性化的一元非线性回归曲线回归
例2 出钢时所用的盛钢水的钢包，由于钢水对耐火材料的侵蚀，
容积不断增大.我们希望知道使用次数与增大的容积之间的关
系.对一钢包作试验，测得的数据列于下表：
使用次数
增大容积
使用次数
增大容积
2
6.42
10
10.49
3
8.20
11
10.59
4
9.58
12
10.60
5
9.50
13
10.80
6
9.70
14
10.60
7
10.00
15
10.90
8
9.93
16
10.76
9
9.99
此即非线性回归或曲线回归问题（需要配曲线）
配曲线的一般方法是：
先对两个变量 x 和 y 作 n 次试验观察得 ( xi , yi ), i  1,2,...,n 画出散
点图，根据散点图确定须配曲线的类型.然后通过变量代换把非线性回归
化成线性回归，即采用非线性回归线性化的方法.
11
10.5
10
散
点
图
9.5
9
8.5
8
7.5
7
6.5
6
2
4
6
8
10
12
14
16
通常选择的六类曲线如下：
1
b
（1）双曲线  a 
y
x
（2）幂函数曲线 y=a x b , 其中 x>0,a>0
（3）指数曲线 y=ae bx 其中参数 a>0.
（4）倒指数曲线 y=ae b / x 其中 a>0，
（5）对数曲线 y=a+blogx,x>0
（6）S 型曲线 y 
解例 2.由散点图我们选配倒指数曲线 y=a e b / x
1
a  be x 根据线性化方法，算得 bˆ  1.1107, Aˆ  2.4587
由此
最后得
ˆ
aˆ  e A  11.6789
y  11.6789e
1.1107

x
二. 多元数据处理方法
1、二维插值
2、多元回归分析
3、聚类分析
4、主成分分析
二维插值的定义
第一种（网格节点）：
y
O















x
已知 mn个节点
其中
互不相同，不妨设
构造一个二元函数
再用
计算插值，即
通过全部已知节点,即
第二种（散乱节点）：
y








0







x
已知n个节点
其中
互不相同，
构造一个二元函数
再用
计算插值，即
通过全部已知节点,即
最邻近插值
y

(x1, y2)




(x1, y1) (x2, y1)

O

(x2, y2)








x
二维或高维情形的最邻近插值，与被插值点最邻近的
节点的函数值即为所求。
注意：最邻近插值一般不连续。具有连续性的最简单
的插值是分片线性插值。
分片线性插值
y

(xi, yj+1) (xi+1, yj+1)







(xi, yj) (xi+1, yj)





O


x
将四个插值点（矩形的四个顶点）处的函数值
依次简记为：
f (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f4
分两片的函数表达式如下：
第一片(下三角形区域)： (x, y)满足
插值函数为：
f ( x, y)  f1  (f 2  f1 )( x  x i )  (f 3  f 2 )( y  y j )
第二片(上三角形区域)：(x, y)满足
y j1  y j
y
(x  x i )  y i
x i 1  x i
插值函数为：
f ( x, y)  f1  (f 4  f1 )( y  y j )  (f 3  f 4 )( x  x i )
注意：(x, y)当然应该是在插值节点所形成的矩形区
域内。显然，分片线性插值函数是连续的；
双线性插值
y

(x1, y2)




(x1, y1) (x2, y1)

O

(x2, y2)








x
双线性插值是一片一片的空间二次曲面构成。
双线性插值函数的形式如下：
f ( x, y)  (ax  b)(cy  d)
其中有四个待定系数，利用该函数在矩形的四个顶
点（插值节点）的函数值，得到四个代数方程，正
好确定四个系数。
用MATLAB作网格节点数据的插值
z=interp2(x0,y0,z0,x,y,’method’)
被插值点
的函数值
插值
节点
被插值点
插值方法
‘nearest’
最邻近插值
‘linear’
双线性插值
‘cubic’
双三次插值
缺省时, 双线性插值
要求x0,y0单调；x，y可取为矩阵，或x取
行向量，y取为列向量，x,y的值分别不能超出
x0,y0的范围。
例：测得平板表面3*5网格点处的温度分别为：
82 81 80 82 84
79 63 61 65 81
84 84 82 85 86
试作出平板表面的温度分布曲面z=f(x,y)的图形。
1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.
输入以下命令：
x=1:5;
y=1:3;
temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];
mesh(x,y,temps)
2．以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.
再输入以下命令:
xi=1:0.2:5;
yi=1:0.2:3;
zi=interp2(x,y,temps,xi',yi,'cubic');
mesh(xi,yi,zi)
画出插值后的温度分布曲面图.
例 山区地貌：
在某山区测得一些地点的高程如下表。平面区域为
1200<=x<=4000,1200<=y<=3600)
试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。
X
Y
1200
1600
2000
2400
2800
3200
3600
1200
1600
2000
2400
2800
3200
3600
4000
1130
1320
1390
1500
1500
1500
1480
1250
1450
1500
1200
1200
1550
1500
1280
1420
1500
1100
1100
1600
1550
1230
1400
1400
1350
1550
1550
1510
1040
1300
900
1450
1600
1600
1430
900
700
1100
1200
1550
1600
1300
500
900
1060
1150
1380
1600
1200
700
850
950
1010
1070
1550
980
通过此例对最近邻点插值、双线性插值方法和双三次插值
方法的插值效果进行比较。
用MATLAB作散点数据的插值计算
插值函数griddata格式为:
cz =griddata（x，y，z，cx，cy，‘method’）
被插值点
的函数值
插值
节点
被插值点
插值方法
‘nearest’ 最邻近插值
‘linear’ 双线性插值
‘cubic’
双三次插值
'v4'- Matlab提供的插值方法
缺省时, 双线性插值
要求cx取行向量，cy取为列向量。
例 在某海域测得一些点(x,y)处的水深z由下
表给出，船的吃水深度为5英尺，在矩形区域（75，
200）*（-50，150）里的哪些地方船要避免进入。
x
y
z
129 140 103.5 88 185.5 195
7.5 141.5 23
147 22.5 137.5
4
8
6
8
6
8
x
y
z
157.5
-6.5
9
107.5
-81
9
77
3
8
105
85.5
8
81
162 162 117.5
56.5 -66.5 84 -33.5
8
9
4
9
1.输入插值基点数据.
2.在矩形区域75,200  50,150作二维插值.
三次插值法(hd1)
3、作海底曲面图
4.作出水深小于5的海域范围,即z=5的等高线.
实验作业
山区地貌：在某山区测得一些地点的高程如下表：(平
面区域1200<=x<=4000,1200<=y<=3600)，试作出该山区的
地貌图和等高线图，并对几种插值方法进行比较。
3600
1480 1500 1550 1510 1430 1300 1200 980
3200
1500 1550 1600 1550 1600 1600 1600 1550
2800
1500 1200 1100 1550 1600 1550 1380 1070
2400
1500 1200 1100 1350 1450 1200 1150 1010
2000
1390 1500 1500 1400 900 1100 1060 950
1600
1320 1450 1420 1400 1300 700 900 850
1200
1130 1250 1280 1230 1040 900 500 700
Y/x
1200 1600 2000 2400 2800 3200 3600 4000
1.4 多元线性回归
多元线性回归在工程上更为有用。
1.4.1数学模型及定义
一般地，影响试验指标的因素不只一个，假设
它们之间有如下的线性关系：
y  0  1 x1  ...   k xk＋
(1)
其中y为可观测的随机变量，x1 , x2 , ..., xk为非随机的
可精确观测的变量， 0，1，. . . . ,  k为k  1个未知
参数， 为随机变量，设E ( )  0, D( )   2  0.
为了估计未知参数 0，1，. . . . ,  k和 2，我
们对x1 , x2 , ..., xk 和y作n次观测得n组观测值
（xi1 , xi 2 , ..., xik , yi ）( i =1, 2, 3, . . . , n) .
它们满足关系式：
yi   0  1 xi1  ...   k xik＋ i
( 2)
i =1, 2, 3, . . . n
其中 1 ,  2 , ...,  n相互独立且是与同分布的随机变量。
为了用矩阵表示上式，令
 y1 
 1 x11
 ... 
1 x
21
Y   ，X 
... ...
 ... 

 
y
 1 x n1
 n
x12
x 22
...
xn2
... x1k 
 0 
 1 
 
 
... x 2 k 
 ，    1  ，   2 
 ... 
 ...
... ... 

 
 
... x nk 

 k
 n 
于是，（）式变为：
1
Y  X 
(3)
其中X 为n  (k  1)阶矩阵，称为资料矩阵， 为k  1维
未知的列向量， 满足：E  0，COV ( ,  )   2 I n .
Y  X  

2
E
(

)

0
,
COV
(

,

)


In

一般称
为 k 元线性回归模型，并简记为 (Y , X , 2 I n )
对（）式取期望，有
1
y  0  1x1  ...  k xk
称为回归平面方程。
线性模型 (Y , X ,  2 I n ) 考虑的主要问题是：
（1） 用试验值（样本值）对未知参数  和  2 作点估计和假设检验，
从而建立 y 与 x1 , x2 ,...,xk 之间的数量关系；
(2) 在 x1  x01 , x2  x02 ,...,xk  x0k , 处对 y 的值作预测与控制，并对 y
作区间估计.
1.4.2模型参数估计
1.对i和 2作估计
用最小二乘法求  0 ,...,  k 的估计量：作离差平方和
n
Q    yi   0  1 xi1  ...   k xik 
i 1
选择  0 ,...,  k 使 Q 达到最小。
根据微积分求极值的方法，得正规方程组
2
 Q
   0
0

 Q
0

  0
............

 Q
   0
k

解得估计值

T
ˆ
 X X
 X Y 
1
T
得到的 ˆ i 代入回归平面方程得：
y  ˆ0  ˆ1 x1  ...  ˆ k xk
称为经验回归平面方程. ˆ i 称为经验回归系
数.
注意：ˆ服从k  1维正态分布，且为的无偏估计，
协方差阵为 2 ( X T X )1.
1.4.3 多项式回归
设变量 x、Y 的回归模型为
Y   m xm  ...   2 x2  1x  0  
其中 m 是已知的，  i (i  1,2,, p) 是未知参数，  服从正
态分布 N (0, 2 ) .
Y  m xm  ...   2 x2  1x  0
称为回归多项式.上面的回归模型称为多项式回归.
令 xi  xi ，i=1，2，…，k 多项式回归模型变
为多元线性回归模型.
1.4.4多元线性回归中的检验与预测
（1）线性模型和回归系数的检验
现在我们来考虑如下模型检验问题
Y  X   

2
 N n (0,  I n )
在实际问题中，y   0  1 x1  ...   k xk 往往只是
一种假设，在求出线性回归方程之后，还须对求出的
线性回归方程同实际观测数据拟合效果进行检验。
可提出以下原假设
H:
0 1  ....  k  0
F检验法
当 H0 成立时， F 
U /k
Qe /( n  k  1)
~ F ( k , n  k  1)
如果 F > F1-α（k，n-k-1）
，则拒绝 H0 ，认为 y 与 x1 ,…, xk 之间显著
地有线性关系；否则就接受
H0 ，认为 y 与 x1,…,
xk 之间线性关系不
n
n
显著. 其中 U 
 yˆ i  y 2 （回归平方和） Qe   ( yi  yˆ i ) 2 (残差平方和）

i 1
i 1
(2)预测
（A）点预测
求出回归方程 yˆ  ˆ0  ˆ1 x1  ... ˆk xk ，对于给定自
*
*
*

*
ˆ
ˆ
ˆ
ˆ
变量的值 x ,...,x ，用 y     x  ...  x 来预测
1
k
0
1 1
k
k
*
*
y    0  1 x1  ...  k xk   .称 yˆ * 为 y * 的点预测.
（B）区间预测
y 的 1   的预测区间（置信）区间为 ( yˆ1 , yˆ 2 ) ,其中
k
k

 yˆ1  yˆ  ˆ e 1   cij xi x j t1 / 2 (n  k  1)
i 0 j 0


k
k
ˆ
 y2  yˆ  ˆ e 1   cij xi x j t1 / 2 (n  k  1)
i 0 j 0

-1
C=L =(cij), L=X’X
Qe
ˆ e 
n  k 1
1.5 逐步回归分析
实际问题中影响因变量的因素可能很多，我们希望从中
挑选出影响显著的自变量来建立回归模型，这就涉及到变量
选择的问题。逐步回归是一种从众多变量中有效地选择重要
变量的方法。它是在多元线性回归的基础上派生出来的一种
算法技巧。
如果采用的自变量越多，则回归平方和越大，残差平方
和越小，然而较多的变量来拟合回归方程，得到的防策划能
够稳定性差，用它作预测可靠性差，精度低．另一方面，如
果采用了y 影响较小的变量而遗漏了重要变量，可导致估计量
产生偏崎和不一致性．为此，我们希望得到“最优”的回归
方程．
“最优”的回归方程就是包含所有对Y有影响的变量, 而不
包含对Y影响不显著的变量回归方程。
选择“最优”的回归方程有以下几种方法：
（1）从所有可能的因子（变量）组合的回归方程中选择最优者；
（2）从包含全部变量的回归方程中逐次剔除不显著因子；
（3）从一个变量开始，把变量逐个引入方程；
（4）“有进有出”的逐步回归分析。
以第四种方法，即逐步回归分析法在筛选变量方面较为理想.
逐步回归分析法的思想：
• 从一个自变量开始，视自变量Y作用的显著程度，从大
到小地依次逐个引入回归方程。
• 当引入的自变量由于后面变量的引入而变得不显著时，要
将其剔除掉。
• 引入一个自变量或从回归方程中剔除一个自变量，为
逐步回归的一步。
• 对于每一步都要进行Y值检验，以确保每次引入新的显
著性变量前回归方程中只包含对Y作用显著的变量。
• 这个过程反复进行，直至既无不显著的变量从回归方程中
剔除，又无显著变量可引入回归方程时为止。
1.６ MATLAB统计工具箱中的回归分析命令
1.６.1多元线性回归
y   0  1 x1  ...   p x p
1)确定回归系数的点估计值：
b=regress( Y, X )
 ˆ 0 
 ˆ 

b   1 
...


ˆ

 p 
 Y1 
Y 
Y   2
 ... 
 
Yn 
对一元线性回归，取p=1即可.
 1 x11
1 x
21
X 
... ...

 1 x n1
x12
x 22
...
xn 2
... x1 p 
... x 2 p 
... ... 

... x np 
2)求回归系数的点估计和区间估计、并检验回归模型：
[b, bint,r,rint,stats]=regress(Y,X,alpha)
残
差
置
信
区
间
用于检验回归模型的统计量，
有三个数值：
相关系数r2、F值、与F对应的
概率p
（
缺
省显
时著
为性
0 水
平
05
）
.
回
归
系
数
的
区
间
估
计
相关系数 r2 越接近 1，说明回归方程越显著；
F > F1-α（k，n-k-1）时拒绝 H0，F 越大，说明回归方程越显著；
与 F 对应的概率 p  时拒绝 H0 ，回归模型成立.
3、画出残差及其置信区间：
rcoplot（r，rint）
例1 解：1、输入数据：
x=[143 145 146 147 149 150 153 154 155 156 157 158 159
160 162 164]';
题目
X=[ones(16,1) x];
Y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]';
2、回归分析及检验：
[b,bint,r,rint,stats]=regress(Y,X)
b,bint,stats
得结果：b =
bint =
-16.0730
-33.7071
0.7194
0.6047
To MATLAB(liti11)
1.5612
0.8340
stats =
0.9282 180.9531
0.0000
即 ˆ 0  16.073, ˆ1  0.7194； ˆ 0 的置信区间为[-33.7017，1.5612], ˆ1 的置信区间为[0.6047,0.834];
r2=0.9282, F=180.9531, p=0.0000
p<0.05, 可知回归模型 y=-16.073+0.7194x 成立.
3、残差分析，作残差图：
rcoplot(r,rint)
从残差图可以看出，除第二个数据外，其余数据的残
差离零点均较近，且残差的置信区间均包含零点，这说明
回归模型 y=-16.073+0.7194x能较好的符合原始数据，而第
二个数据可视为异常点.
Residual Case Order Plot
4
2
1
Residuals
4、预测及作图：
z=b(1)+b(2)*x
plot(x,Y,'k+',x,z,'r')
3
0
-1
-2
-3
-4
-5
2
4
6
8
10
Case Number
12
14
16
To MATLAB(liti12)
返回
1.６.2多 项 式 回 归
（1）一元多项式回归
A、回归：
y=a1xm+a2xm-1+…+amx+am+1
1）确定多项式系数的命令：[p,S]=polyfit（x，y，m）
其中 x=（x1，x2，…，xn）
，y=（y1，y2，…，yn）
；p=
（a1，a2，…，am+1）是多项式 y=a1xm+a2xm-1+…+amx+am+1
的系数；S 是一个矩阵，用来估计预测误差.用于其它
函数的计算
2）一元多项式回归命令：polytool（x，y，m）
此命令产生一个交互式的画面，画面中有拟合曲线和y
的置信区间。通过左下方的Export菜单，可以输出回归系数
等。
B、预测和预测误差估计：
（1）Y=polyval（p，x）求polyfit所得的回归多项式在x处的预
测值Y；
（2）[Y，DELTA]=polyconf（p，x，S，alpha）求polyfit所得
的回归多项式在x处的预测值Y及预测值的显著性为1alpha的置信区间Y DELTA；alpha缺省时为0.05.
一元多项式回归也可以化为多元线性回归来解。
例 2 观测物体降落的距离 s 与时间 t 的关系，得到数据如下表，求 s
关于 t 的回归方程sˆ  a  bt  ct 2 .
t
s
(s)
1/30
2/30
3/30
4/30
5/30
6/30
7/30
(cm)
11.86
15.67
20.60
26.69
33.71
41.93
51.13
(s)
8/30
9/30
10/30
11/30
12/30
13/30
14/30
(cm)
61.49
72.90
85.44
99.08
113.77
129.54
146.48
t
s
法一
直接作二次多项式回归：
t=1/30:1/30:14/30;
s=[11.86 15.67 20.60 26.69 33.71 41.93 51.13 61.49 72.90
85.44 99.08 113.77 129.54 146.48];
[p,S]=polyfit(t,s,2)
To MATLAB（liti21）
得回归模型为 ：
sˆ  489.2946t 2  65.8896t  9.1329
法二
化为多元线性回归：
To MATLAB(liti22)
t=1/30:1/30:14/30;
s=[11.86 15.67 20.60 26.69 33.71 41.93 51.13 61.49 72.90
85.44 99.08 113.77 129.54 146.48];
T=[ones(14,1) t' (t.^2)'];
[b,bint,r,rint,stats]=regress(s',T);
b,stats
得回归模型为 ：
sˆ  9.1329 65.8896t  489.2946t 2
预测及作图
Y=polyconf(p,t,S)
plot(t,s,'k+',t,Y,'r')
To MATLAB(liti23)
1.6.3多元二项式回归
命令：rstool（x，y，’model’, alpha）
nm矩阵
显著性水平
（缺省时为0.05）
n维列向量
由下列 4 个模型中选择 1 个（用字符串输入，缺省时为线性模型）：
linear（线性）： y  0  1x1 
 mxm
purequadratic（纯二次）： y   0  1 x1 
interaction（交叉）： y   0  1 x1 
  m x m    jjx 2j
 mxm 
quadratic（完全二次）： y   0  1x1 
n

1 j k  m
 mxm 
j1
 jk x jx k

1 j,k  m
 jk x jx k
命令rstool产生一个交互式画面，画面中有m个图形，这m
个图形分别给出了一个独立变量xi（另m-1个变量取固定值）
与y的拟合曲线，以及y的置信区间。可以通过键入不同的xi值
来获得相应的y值。
例3 设某商品的需求量与消费者的平均收入、商品价格的统计数
据如下，建立回归模型，预测平均收入为800、价格为6时
的商品需求量.
需求量
收入
价格
100
1000
5
75
600
7
80
1200
6
70
500
6
50
300
8
65
400
7
90
1300
5
100
1100
4
110
1300
3
60
300
9
选择纯二次模型，即 y   0  1 x1   2 x 2  11 x12   22 x 22
解
直接用多元二项式回归：
x1=[1000 600 1200 500 300 400 1300 1100 1300 300];
x2=[5 7 6 6 8 7 5 4 3 9];
y=[100 75 80 70 50 65 90 100 110 60]';
x=[x1' x2'];
rstool(x,y,'purequadratic')
在左边图形下方的方框中输入800，右边图形下方的方框中输入6。
则画面左边的“Predicted Y”下方的数据变为86.3971，即预测出平均收
入为800、价格为6时的商品需求量为86.3971.
在画面左下方的下拉式菜单中选”all”, 则beta（回归系数）、
rmse（剩余标准差）和residuals（残差）都传送到Matlab工作区中.
在Matlab工作区中输入命令： beta, rmse
得结果：beta =
110.5313
0.1464
-26.5709
-0.0001
1.8475
rmse =
4.5362
故回归模型为： y  110.5313 0.1464x1  26.5709x 2  0.0001x12  1.8475x 22
剩余标准差为 4.5362, 说明此回归模型的显著性较好.
To MATLAB(liti31)
1.７非线性回 归
1.７.1回归：
是事先用m-文件定
义的非线性函数
（1）确定回归系数的命令：
[beta，r，J]=nlinfit（x，y，’model’, beta0）
估计出的
回归系数
残差
Jacobian矩阵，用于
估计预测误差需要
的数据。
输入数据x、y分别为
n m矩阵和n维列向
量，对一元非线性回
归，x为n维列向量。
回归系数
的初值
（2）非线性回归命令：nlintool（x，y，’model’, beta0，alpha）
其中个参数含义同前，alpha为显著性水平，缺省时为0.05。
该命令产生一个交互式的画面，画面中有拟合曲线和y的置信
区间。通过左下方的Export菜单，可以输出回归系数等。
1.７.2、预测和预测误差估计：
[Y，DELTA]=nlpredci（’model’, x，beta，r，J）
该命令用于求nlinfit 或nlintool所得的回归函数在x处的预
测值Y及预测值的显著性为1-alpha的置信区间Y DELTA.
例 4 对第一节例2，求解如下：
题目
1、对将要拟合的非线性模型 y=ae b / x ，建立 m-文件 volum.m 如下：
function yhat=volum(beta,x)
yhat=beta(1)*exp(beta(2)./x);
2、输入数据：
x=2:16;
y=[6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60
10.90 10.76];
beta0=[8 2]';
To MATLAB(liti41)
3、求回归系数：
[beta,r ,J]=nlinfit(x',y','volum',beta0)；
beta
即得回归模型为：
得结果：beta =
1.10641
11.6036

x
y

11
.
6036
e
-1.0641
4、预测及作图：
[YY,delta]=nlpredci('volum',x',beta,r ,J)；
plot(x,y,'k+',x,YY,'r')
１.８逐 步 回 归
显著性水平（缺省时为0.０5）
逐步回归的命令是：
stepwise（x，y，inmodel，alpha）
自变量数据,
n m阶矩阵
因变量数据，
n 1 阶矩阵
矩阵的列数的指标，给出初
始模型中包括的子集（缺省
时设定为全部自变量）
运行stepwise命令时产生三个图形窗口：Stepwise Plot，
Stepwise Table，Stepwise History.
在Stepwise Plot窗口，显示出各项的回归系数及其置信区间.
Stepwise Table 窗口中列出了一个统计表，包括回归系数及
其置信区间，以及模型的统计量剩余标准差（RMSE）、相关系
数（R-square）、F值、与F对应的概率P.
例6 水泥凝固时放出的热量y与水泥中4种化学成分x1、x2、x3、 x4
有关，今测得一组数据如下，试用逐步回归法确定一个 线性模
型.
序号
1
2
3
4
5
6
7
8
9
10
11
12
13
x1
x2
x3
x4
7
26
6
60
1
29
15
52
11
56
8
20
11
31
8
47
7
52
6
33
11
55
9
22
3
71
17
6
1
31
22
44
2
54
18
22
21
47
4
26
1
40
23
34
11
66
9
12
10
68
8
12
y
78.5
74.3
104.3
87.6
95.9
109.2
102.7
72.5
93.1
115.9
83.8
113.3
109.4
1、数据输入：
x1=[7 1 11 11 7 11 3 1 2 21 1 11 10]';
x2=[26 29 56 31 52 55 71 31 54 47 40 66 68]';
x3=[6 15 8 8 6 9 17 22 18 4 23 9 8]';
x4=[60 52 20 47 33 22 6 44 22 26 34 12 12]';
y=[78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3
109.4]';
x=[x1 x2 x3 x4];
2、逐步回归：
（1）先在初始模型中取全部自变量：
stepwise(x,y)
得图Stepwise Plot 和表Stepwise Table
图Stepwise Plot中四条直线都是虚
线，说明模型的显著性不好
Stepwise Table中列出了各变量的回
归系数和置信区间，从该表中看出
变量x3和x4的显著性最差.
（2）在图Stepwise Plot中点击直线3和直线4，移去变量x3和x4
移去变量x3和x4后模型具有显著性.
虽然剩余标准差（RMSE）没
有太大的变化，但是统计量F的
值明显增大，因此新的回归模型
更好.
To MATLAB(liti51）
（3）对变量y和x1、x2作线性回归：
X=[ones(13,1) x1 x2];
b=regress(y,X)
得结果：b =
52.5773
1.4683
0.6623
故最终模型为：y=52.5773+1.4683x1+0.6623x2
3、聚类分析
§1 什么是聚类分析
例 对10位应聘者做智能检验。3项指标X，Y
和Z分别表示数学推理能力，空间想象能力和语
言理解能力。其得分如下，选择合适的统计方
法对应聘者进行分类。
应聘者
X
1
28
2
18
3
11
4
21
5
26
6
20
7
16
8
14
9
24
10
22
Y
29
23
22
23
29
23
22
23
29
27
Z
28
18
16
22
26
22
22
24
24
24
应聘者
X
1
28
2 3 4 5 6 7 8 9 10
18 11 21 26 20 16 14 24 22
Y
Z
29
28
23 22 23 29 23 22 23 29 27
18 16 22 26 22 22 24 24 24
我们直观地来看，这个分类是否合理？
计算4号和6号得分的离差平方和：
(21-20)2+(23-23)2+(22-22)2=1
计算1号和2号得分的离差平方和：
(28-18)2+(29-23)2+(28-18)2=236
计算1号和3号得分的离差平方和为482，由此可
见一般，分类可能是合理的，欧氏距离很大的应聘
者没有被聚在一起。
由此，我们的问题是如何来选择样品间相似的
测度指标，如何将有相似性的类连接起来？
聚类分析根据一批样品的许多观测指标，按照
一定的数学公式具体地计算一些样品或一些参数
(指标)的相似程度，把相似的样品或指标归为一类，
把不相似的归为一类。
例如对上市公司的经营业绩进行分类；据经
济信息和市场行情，客观地对不同商品、不同用
户及时地进行分类。又例如当我们对企业的经济
效益进行评价时，建立了一个由多个指标组成的
指标体系，由于信息的重叠，一些指标之间存在
很强的相关性，所以需要将相似的指标聚为一类，
从而达到简化指标体系的目的。
§2 相似系数和距离
一、变量测量尺度的类型
为了将样本进行分类，就需要研究样品之间的关
系；而为了将变量进行分类，就需要研究变量之间
的关系。但无论是样品之间的关系，还是变量之间
的关系，都是用变量来描述的，变量的类型不同，
描述方法也就不同。通常，变量按照测量它们的尺
度不同，可以分为三类。
(1)间隔尺度。指标度量时用数量来表示，其数
值由测量或计数、统计得到，如长度、重量、收入、
支出等。一般来说，计数得到的数量是离散数量，
测量得到的数量是连续数量。在间隔尺度中如果存
在绝对零点，又称比例尺度。
(2)顺序尺度。指标度量时没有明确的数量表示，只
有次序关系，或虽用数量表示，但相邻两数值之间的差距
并不相等，它只表示一个有序状态序列。如评价酒的味道，
分成好、中、次三等，三等有次序关系，但没有数量表示。
又如评价产品的质量，虽可用一、二、三等来表示，但一
等与二等之间、二等与三等之间的差距并不一定相等。
(3)名义尺度。指标度量时既没有数量表示也没有次
序关系，只有一些特性状态，如眼睛的颜色，化学中催化
剂的种类等。在名义尺度中只取两种特性状态的变量是很
重要的，如电路的开和关，天气的有雨和无雨，人口性别
的男和女，医疗诊断中的“十”和“一”，市场交易中的
买和卖等都是此类变量。显然，对于具有多个特性状态的
变量，可通过并类的方法将其转化为二性状态变量。
二、数据的变换处理
所谓数据变换，就是将原始数据矩阵中的每个元素，
按照某种特定的运算把它变成为一个新值，而且数值的变
化不依赖于原始数据集合中其它数据的新值。
1、中心化变换
中心化变换是一种坐标轴平移处理方法，它是先求出
每个变量的样本平均值，再从原始数据中减去该变量的均
值，就得到中心化变换后的数据。
设原始观测数据矩阵为：
 x11 x12  x1 p 
x

x

x
21
22
2p

X

 

x

x

x
np 
 n1 n 2
xij*  xij  x j
(i  1,2,3,, n; j  1,2,3,, p)
中心化变换的结果是使每列数据之和
均为0，即每个变量的均值为0，而且每列
数据的平方和是该列变量样本方差的(n—1)
倍，任何不同两列数据之交叉乘积是这两
列变量样本协方差的(n—1)倍，所以这是
一种很方便地计算方差与协方差的变换。
2、极差规格化变换
规格化变换是从数据矩阵的每一个变量中找出其
最大值和最小值，这两者之差称为极差，然后从每
个变量的每个原始数据中减去该变量中的最小值，
再除以极差，就得到规格化数据。即有：
xij* 
xij  min( xij )
i 1, 2 ,, n
(i  1,2,3,, n; j  1,2,3,, p)
Rj
R j  max ( xij )  min( xij )
i 1, 2 ,,n
0  xij*  1
i 1, 2 ,,n
经过规格化变换后，数据矩阵中每列即每个变量的最
大数值为1，最小数值为0，其余数据取值均在0－1之间；
并且变换后的数据都不再具有量纲，便于不同的变量之间
的比较。
3、标准化变换
标准化变换也是对变量的数值和量纲进行类似于规格
化变换的一种数据处理方法。首先对每个变量进行中心化
变换，然后用该变量的标准差进行标准化。即有：
xij  x j
(i  1,2,3,, n; j  1,2,3,, p)
x 
Sj
1 n
2
 ( xij  x j )
Sj 
n  1 i 1
*
ij
应聘者
X
1
28
2 3 4 5 6 7 8 9 10
18 11 21 26 20 16 14 24 22
Y
Z
29
28
23 22 23 29 23 22 23 29 27
18 16 22 26 22 22 24 24 24
clear
sj=[28 18 11 21 26 20 16 14 24 22;
29,23,22 23,29 23,22,23,29,27;
28,18,16,22,26,22,22,24,24,24]’;
[a,b]=size(sj);
for j=1:b
jz(j)= mean(sj(:,j))
S(j)=sum((sj(:,j)-jz(j)).^2)/(a-1);
sj1(:,j)=(sj(:,j)-jz(j))/S(j)
end
sj1
1 n
2

Sj 
( xij  x j )
i

1
n 1
xij  x j
x 
Sj
*
ij
经过标准化变换处理后，每个变量即数据矩阵
中每列数据的平均值为0，方差为1，且也不再具有
量纲，同样也便于不同变量之间的比较。变换后，
数据短阵中任何两列数据乘积之和是两个变量相关
系数的（n－1）倍，所以这是一种很方便地计算相
关矩阵的变换。
4．对数变换
对数变换是将各个原始数据取对数，将原始数据的
对数值作为变换后的新值。即：
xij*  log( xij )
三、样品间亲疏程度的测度
研究样品或变量的亲疏程度的数量指标有两
种，一种叫相似系数，性质越接近的变量或样
品，它们的相似系数越接近于1或-l，而彼此
无关的变量或样品它们的相似系数则越接近于
0，相似的为一类，不相似的为不同类；另一
种叫距离，它是将每一个样品看作p维空间的
一个点，并用某种度量测量点与点之间的距离，
距离较近的归为一类，距离较远的点应属于不
同的类。变量之间的聚类即只型聚类分析，常
用相似系数来测度变量之间的亲疏程度。而样
品之间的聚类即Q型聚类分析，则常用距离来
测度样品之间的亲疏程度。
1、定义距离的准则
定义距离要求满足第i个和第j个样品之间的距
离如下四个条件:
dij  0对一切的i和j成立;
dij  0当且仅当i  j成立;
dij  d ji 0对一切的i和j成立;
dij  dik  d kj对于一切的i和j成立.
2、常用距离的算法
(1) 明氏距离测度
设 x i  xi1 , xi 2 , , xip  和 x j  ( x j1 , x j 2 ,, x jp )
是第i和 j 个样品的观测值，则二者之间的距离
为：
明氏距离
特别，欧氏距离
dij  ( k1| xik  x jk |g )
p
dij 
p
 ( xik
k 1
 x jk )
2
1
g
该距离的缺点明考夫斯基距离主要有以下两
个缺点：
①明氏距离的值与各指标的量纲有关，而各
指标计量单位的选择有一定的人为性和随意性，
各变量计量单位的不同不仅使此距离的实际意
义难以说清，而且，任何一个变量计量单位的
改变都会使此距离的数值改变从而使该距离的
数值依赖于各变量计量单位的选择。
②明氏距离的定义没有考虑各个变量之间的
相关性和重要性。实际上，明考夫斯基距离是
把各个变量都同等看待，将两个样品在各个变
量上的离差简单地进行了综合。
(2)杰氏距离
这是杰斐瑞和马突斯塔(Jffreys 8L Matusita)所定义的一种
距离，其计算公式为：
d ij ( J ) 

p

k 1
( xik  x jk )
2

12
(3)兰氏距离
这是兰思和维廉姆斯(Lance SL Williams)所给定的一种距
离，其计算公式为：
p xik  x jk
dij ( L)  
k 1 x  x
ik
jk
这是一个自身标准化的量，由于它对大的奇异值不敏感
，这样使得它特别适合于高度偏倚的数据。虽然这个距离
有助于克服明氏距离的第一个缺点，但它也没有考虑指标
之间的相关性。
(4)马氏距离
这是印度著名统计学家马哈拉诺比斯
(P．C．’Mahalanobis)所定义的一种距离，
其计算公式为：
dij2  (xi  x j ) 1 (xi  x j )
分别表示第i个样品和第j样品的p指标观测值所组成的列
向量，即样本数据矩阵中第i个和第j个行向量的转置，表
示观测变量之间的协方差短阵。在实践应用中，若总体协
方差矩阵未知，则可用样本协方差矩阵作为估计代替计
算。
应聘者
1
2
3
4
5
6
7
8
9
10
X
28
18
11
21
26
20
16
14
24
22
Y
29
23
22
23
29
23
22
23
29
27
Z
28
18
16
22
26
22
22
24
24
24
dij2  (xi  x j ) 1 (xi  x j )
标准化
X
Y
Z
1
0.2791
0.4186
0.4324
2
-0.0698
-0.2093
-0.3683
for j=1:a
3
-0.3140
-0.3140
-0.5285
4
0.0349
-0.2093
-0.0480
D(i,j)=(sj(i,:)-sj(j,:))*var1* (sj(i,:)-sj(j,:))’;
5
0.2093
0.4186
0.2722
6
0
-0.2093
-0.0480
7
-0.1395
-0.3140
-0.0480
end
8
-0.2093
-0.2093
0.1121
9
0.1395
0.4186
0.1121
D
10
0.0698
0.2093
0.1121
var1=cov(sj1);
for i=1:a
end
马氏距离又称为广义欧氏距离。显然，马氏距离与上述
各种距离的主要不同就是马氏距离考虑了观测变量之间的
相关性。如果假定各变量之间相互独立，即观测变量的协
方差矩阵是对角矩阵，则马氏距离就退化为用各个观测指
标的标准差的倒数作为权数进行加权的欧氏距离。因此，
马氏距离不仅考虑了观测变量之间的相关性，而且也考虑
到了各个观测指标取值的差异程度，消除了各个观测指标
不同量纲的影响。这表明，马氏距离对任何非奇异线性变
换都具有不变性。为了对马氏距离和欧氏距离进行一下比
较，以便更清楚地看清二者的区别和联系，现考虑一个例
子。
例如，假设有一个二维正态总体，它的分布为：
 0   1 0.9 
N 2  , 

 0   0.9 1 
 0.9 
1  1
 


0.19   0.9
1 
1
设A(1,1)和B  (1,1)两点。
d A (M )  1.05
d A (U )  2
d B (M )  20
d B (U )  2
马氏距离虽然考虑了观测变量之间的相关性，并且
也不受观测变量量纲不同的影响，但是在聚类分析之
前，如果用全部数据计算的均值向量和协方差矩阵来
计算马氏距离，效果并不是很好。比较合理的办法是
用各个类的样本来计算各自的协方差矩阵，同一类样
本的马氏距离应当用这一类的协方差短阵来计算。然
而，类的形成要依赖于样品之间的距离，反过来样品
间合理的马氏距离又依赖于类，这就形成了一个恶性
循环。因此，在实际聚类分析处理中，马氏距离也不
是理想的距离。
。
5、斜交空间距离
由于各变量之间往往存在着不同的相关
关系，用正交空间的距离来计算样本间的
距离易变形，所以可以采用斜交空间距离。
1 p p

dij   2   ( xih  x jh )( xik  x jk ) hk 
 p h1k 1

12
当各变量之间不相关时，斜交空间退化
为欧氏距离。
四、变量间亲疏程度的测度
1、准则
当对变量进行聚类时，最常见的相似系数
是样本相关系数。性质越接近的变量，它们的
相关系数的绝对值越接近1，而彼此无关的变
量的相关系数则越接近0。满足的条件：
Cij  1,当且仅当两变量之间为完全线性关系；
| Cij | 1, 对一切的i和j成立；
Cij  C ji , 对一切的i和j成立；
2、相似系数的算法
（1）相似系数

设 x i  xi1 , xi 2 , , xip  和 x j  ( x j1 , x j 2 ,, x jp )
是第 i 和 j 个样品的观测值，则二者之间的相似
测度为:
p
其中
 ij 
 ( xik
k 1
 xi )( x jk  x j )
[ k1( xik  xi ) ][ k1( x jk  x j ) 2 ]
p
2
p
（2）夹角余弦
夹角余弦时从向量集合的角度所定义的一种测度变量
之间亲疏程度的相似系数。设在n维空间的向量


x j  x1 j , x2 j ,, xnj 
xi   x1i , x2i ,, xni 
n
cij  cos ij 

xki xkj
k 1
n
2 n
 ki 
k 1
k 1
x
xkj2
观测向量在原点出的夹角的余弦，若变量之间的
相关程度密切，则夹角接近0，其余弦接近1；反之
接近0。
dij2  1  Cij2
五、距离和相似系数选择的原则
一般说来，同一批数据采用不同的亲疏测度指标，会得
到不同的分类结果。产生不同结果的原因，主要是由于不
同的亲疏测度指标所衡量的亲疏程度的实际意义不同，也
就是说，不同的亲疏测度指标代表了不同意义上的亲疏程
度。因此我们在进行聚类分析时，应注意亲疏测度指标的
选择。通常，选择亲疏测度指标时，应注意遵循的基本原
则主要有：
(1)所选择的亲疏测度指标在实际应用中应有明确的意义
。如在经济变量分析中，常用相关系数表示经济变量之间
的亲疏程度。
(2)亲疏测度指标的选择要综合考虑已对样本观测数据
实施了的变换方法和将要采用的聚类分析方法。如在标准
化变换之下，夹角余弦实际上就是相关系数；又如若在进
行聚类分析之前已经对变量的相关性作了处理，则通常就
可采用欧氏距离，而不必选用斜交空间距离。此外，所选
择的亲疏测度指标，还须和所选用的聚类分析方法一致。
如聚类方法若选用离差平方和法，则距离只能选 用欧氏
距离。
(3)适当地考虑计算工作量的大小。如对大样
本的聚类问题，不适宜选择斜交空间距离，因采用
该距离处理时，计算工作量太大。样品间或变量间
亲疏测度指标的选择是一个比较复杂且带主规性的
问题，我们应根据研究对象的特点作具体分折，以
选择出合适的亲疏测度指标。实践中，在开始进行
聚类分析时，不妨试探性地多选择几个亲疏测度指
标，分别进行聚类，然后对聚类分析的结果进行对
比分析，以确定出合适的亲疏测度指标。
六、类的连接方法
x11•
单连接
（最短距离） x12•
完全连接
•
（最长距离）x11•
• •
平均连接
（平均距离）
d13
•
•
x21•
x22•
•
d12
• •
x21•
•
•
d1    d 9
9
•
•
§ 系统聚类方法
（一）方法
1、根据样品的特征，规定样品之间的距离 ，
共有
个。将所有列表，记为D（0）表，该表
Cn2
d ij
是一张对称表。
2、选择D（0）表中最小的非零数，不妨假设 d pq ，
于是将 G p 和 Gq 合并为一类，记为

Gr  Gp，G。
q
3、利用递推公式计算新类与其它类之间的距离。
分别删除D（0）表的第p，q行和第p，q列，并新
增一行和一列添上的结果，产生D（1）表。
G1
G2

Gn
G1
0
d12

d1n
G2
d 21
0

d2n





Gn
d n1
dn2
Gp
Gq

0
（二）常用的种类
1、 最短距离法
设抽取五个样品，每个样品只有一个变量，它们
是1，2，3.5，7，9。用最短距离法对5个样品进行分类。
首先采用绝对距离计算距离矩阵：
D(0)
G1
G2
G3
G4
G1
0
G2
G3
G4
G5
1
0
2.5
1.5
0
6
5
3.5
0
8
7
5.5
2
G5
0
然后 G1 和 G2 被聚为新类 G6 ，得 D (1) ：
G6
G6
G3
G4
G5
G3
G4
G5
0
1.5
0
5
3.5
0
7
5.5
2
定义距离：D pq  Mindij：xi  G p，x j  Gq 
递推公式：Drl  MinD pl，Dql 
l  p，q
0
cophenet函数
调用格式：c=cophenetic(Z,Y)
说明：利用pdist函数生成的Y和linkage函数生成的Z计算cophenet相关系数。
pdist函数
调用格式：Y=pdist(X,’metric’)
说明：用 ‘metric’指定的方法计算 X 数据矩阵中对象之间的距离。
X：一个m×n的矩阵，它是由m个对象组成的数据集，每个对象的大小为n。
metric’取值如下：
‘euclidean’：欧氏距离（默认）；‘seuclidean’：标准化欧氏距离；
‘mahalanobis’：马氏距离；‘cityblock’：布洛克距离；
‘minkowski’：明可夫斯基距离；‘cosine’：
‘chebychev’：Chebychev距离。
linkage函数
调用格式：Z=linkage(Y,’method’)
说 明：用‘method’参数指定的算法计算系统聚类树。
Y：pdist函数返回的距离向量；
method：可取值如下：
‘single’：最短距离法（默认）； ‘complete’：最长距离法；
‘average’：未加权平均距离法； ‘weighted’： 加权平均法；
‘centroid’：质心距离法；
‘median’：加权质心距离法；
‘ward’：内平方距离法（最小方差算法）
应聘者
1
2
3
4
5
6
7
8
9
10
X
28
18
11
21
26
20
16
14
24
22
Y
29
23
22
23
29
23
22
23
29
27
Z
28
18
16
22
26
22
22
24
24
24
for j=1:a-1
for i=j+1:a
Y2(h)=D(i,j);%取出所需距离(10*10-10)/2
h=h+1;
end
end
Z2=linkage(Y2);%定义变量之间的连接
C2=cophenet(Z2,Y2); %评价聚类信息
T=cluster(Z2,2);%表示分成两类
H=dendrogram(Z2);%画出聚类图
clear
sj=[28 18 11 21 26 20 16 14 24 22;
29,23,22 23,29 23,22,23,29,27;
28,18,16,22,26,22,22,24,24,24]'
X2=zscore(sj); %标准化数据
Y2=pdist(X2);%计算距离
Z2=linkage(Y2);
C2=cophenet(Z2,Y2);
T=cluster(Z2,4);%分成4类
figure(2)
H=dendrogram(Z2);%画出聚类图
2、最长距离法
用最长距离法对5个样品进行分类。首先采用绝
对距离计算距离矩阵：
D(0)
G1
G2
G3
G4
G1
0
G2
G3
G4
G5
1
0
2.5
1.5
0
6
5
3.5
0
8
7
5.5
2
G5
0
然后和被聚为新类，得：
G6
G3
G4
G6
0
G3
2.5
0
G4
6
3.5
0
G5
8
5.5
2
定义距离：D pq  Max{dij：Xi  G p，X j  Gq }
递推公式：Drl  Max{Dpl，Dql }，l  p，q
G5
0
3、中间距离法
用中间距离法对5个样品进行分类。首先采用绝
D(0)
对距离计算距离平方矩阵：
G1
G2
G3
G4
G1
0
G2
1
0
G3
6.25
2.25
0
G4
G5
36
25
12.25
0
64
49
30.25
4
G5
0
1 2 1 2 1 2
递推公式：D  Dkp  Dkq  D pq
2
2
4
1 2 1 2
1
2
2
递推公式：Dkr  Dkp  Dkq  D pq ，  〈〈0
2
2
4
2
kr
D13  D23 1
D63 
 D12
2
4
6.25  2.25 1
D63 
 1  4
2
4
G6
G6
G3
G3
G4
G5
0
4
0
G4
30.25
12.25
0
G5
56.25
30.25
4
0
4、类平均法
类平均法定义类间的距离是两类间样品的距离的平
均数。
1

D 
n p nq x G
2
pq
i
G1
G2

p
x j Gq
dij2
G3
G4
G1
0
G2
1
0
G3
6.25
2.25
0
G4
G5
36
25
12.25
0
64
49
30.25
4
G5
0
然后和被聚为新类，得 D (1) ：
G6
G3
G4
G6
0
G3
4.25
0
G4
30.25
12.25
0
G5
56.25
30.25
4
2
2
n
D

n
D
递推公式：Drk2  p kp q kq
n p  nq
G5
0
5、可变类平均法
类平均法的递推公式中，没有反映Gp 类和Gq 类的
距离有多大，进一步将其改进，加入D2Pq ，并给定系
数<1，则类平均法的递推公式改为：
2
2
n
D

n
D
q qk
2
Drk2  (1   ) p pk
 D pq
n p  nq
用此递推公式进行聚类就是可变类平均法。
6、离差平方和法
G1 和 G2 被聚为新类，重心为 X 6  (1  2) / 2  1.5
如G1 和G2 为一类，则离差平方和
S12  (1  1.5)2  (2  1.5)2  0.5
如 G1和 G3为一类，则离差平方和
S13  (1  2.25) 2  (2  2.25) 2  3.125
类似于方差分析的想法，如果类分得恰当，同类内的
样品之间的离差平方和应较小，而类间的离差平方和应当
较大。
离差平方和法的思路是，当k固定时，选择使S达到
最小的分类。先让n个样品各自成一类，然后缩小一类，
每缩小一类离差平方和就要增大，选择使S2增加最小的两
类合并，直到所有的样品归为一类为止。离差平方和法
定义类间的平方距离为
G1
G1
G2
G3
G4
G5
G2
G3
G4
G5
0
0.5
0
3.125
1.125
0
18
12.50
6.125
0
32
24.50
15.125
2
0
2
定义距离为离差平方和 的增量：Dpq
 Sr2  S p2  Sq2
其中 S 2 是由Gp 和Gq 合并成的Gr 类的类内离差平方和。
r
可以证明离差平方和的聚类公式为
nk  n p 2 nk  nq 2
nk
2
递推公式：D 
D pk 
Dqk 
D pq
nr  nk
nr  nk
nk  nr
2
rk
7、可变方法
如果让中间距离法的递推公式前两项的系数也依赖于，
则递推公式为：
1  2
2
2
D 
( Dkp  Dkq
)  D pq
，〈1
2
2
kr
用上式作为递推公式的系统聚类法称为可变法。
8、重心法
用重心法对5个样品进行分类。首先采用
绝对距离计算距离平方矩阵：
G1
G2
G3
G4
G1
0
G2
1
0
G3
6.25
2.25
0
G4
G5
36
25
12.25
0
64
49
30.25
4
G5
0
重心法，也称为样品的均值法。设Gp和Gq 为两个
1
Xp 
np
类
np

xi G p
xi 和X q
1

nq
nq

xi Gq
xi
分别为Gp和Gq的重心，类与类之间的距离定义为
两个类重心（类内样品平均值）间的平方距离。
G4和G6的距离为
1 2 1 2 1 2 1
1
1
D  D41  D42  D12   36   25   1  30.25
2
2
4
2
2
4
2
46
设某一步Gp和Gq的重心分别为为和，类内的样品
数分别为和，如果要把Gp 和Gq 合并为Gr 类，则Gr 类的
样品数nr=np+nq ，Gr 类的重心为X p 和 X q的加权算术平
均数：
n p X p  nq X q
Xr 
n p  nq
递推公式为：
n p 2 nq 2 n p nq 2
D  D pk  Dqk  2 D pq
nr
nr
nr
2
kr
（三）确定类的个数
在聚类分析过程中类的个数如何来确定才合适呢？这是一个十
分困难的问题，人们至今仍未找到令人满意的方法。但是这个
问题又是不可回避的。下面我们介绍几种方法。
1、给定阈值——通过观测聚类图，给出一个合适的阈值T。
要求类与类之间的距离不要超过T值。例如我们给定T=0.35，当
聚类时，类间的距离已经超过了0.35，则聚类结束。
总离差平方和的分解（准备知识）
 x11
x
 21
 
x
 n1
x1
x12  x1 p 
x22  x2 p 


 
xn 2  xnp 
x2
xp
总离差平方和 （x11  x1 ) 2    ( xn1  x1 ) 2
  ( x1 p  x p ) 2   ( xnp  x p ) 2
如果着些样品被分成两类
 x11
x
 21
 
x
 n1
xn 2
(1)
1
(1)
2
1
x
x12
x22

1
x
 x1 p 
 x2 p 

 
 xn p 
1
x p(1)
 x11
x
 21
 
x
 n1
x1( 2)
2
x1 p 
x2 p 

 
 xn p 


x12
x22

xn 2
2
2
x p( 2 )
x2( 2)
一组的离差平方和 （x11  x（1 1）)2    ( xn 1  x（1 1）)2
1
  ( x1 p  x p(1) )2  ( xn p  x p(1) )2
1
二组的离差平方和 （x11  x（1 2）) 2    ( xn 1  x（1 2）) 2
2
  ( x1 p  x p( 2) ) 2  ( xn p  x p( 2) ) 2
2
可以证明：
总离差平方和
＝组内离差平方和＋组间离差平方和
令T为总离差平方和
令PG为分为G类的组内离差平方和。
PG
2、统计量 R  1 
T
其中T是数据的总离差平方和， P 是组内离差平方和。
G
R 2 比较大，说明分G个类时类内的离差平方和比较小，
也就是说分G类是合适的。但是，分类越多，每个类
的类内的离差平方和就越小， 也就越大；所以我们
R2
只能取合适的G，使得
足够大，而G本生很小，随
R2
着G的增加， 的增幅不大。比如，假定分4类时，
R2
=0.8；下一次合并分三类时，下降了许多，
=0.32，
R2
R2
则分4 类是合适的。
2
3、伪F统计量的定义为 F  (T  PG ) (G  1)
PG (n  G )
伪F统计量用于评价聚为G类的效果。如果
聚类的效果好，类间的离差平方和相对于
类内的离差平方和大，所以应该取伪F统计
量较大而类数较小的聚类水平。
Pseudo F Statistic
120
110
100
90
80
70
60
50
40
30
20
10
0
1
2
3
4
5
6
7
8
9
10
11
Number of Clusters
12
13
14
15
16
17
18
19
4、伪 t 2 统计量的定义为
BKL
t 
(WK  WL ) ( N K  N L  2)
2
其中 WL和 WK分别是的类内离差平方和，
是将K和L合并为第M类的离差平方和
BKL = WM - WL - WK
为合并导致的类内离差平方和的增量。用它
评价合并第K和L类的效果，伪 t 2 统计量大说
明不应该合并这两类，应该取合并前的水平。
五、 系统聚类法的基本性质
（一） 单调性
在聚类分析过程中，并类距离分别为l k （k=1，2，
3，…）若满足 l1  l2    lk  lk 1   ，则称该聚类
方法具有单调性。可以证明除了重心法和中间距离法之外，
其他的系统聚类法均满足单调性的条件。
（二）空间的浓缩和扩张
1、 定义矩阵的大小
设同阶矩阵D（A）和D（B），如果D（A）的每一个元素
不小于D（B）的每一个元素，则记为 D( A)  D( B)。
2、空间的浓缩和扩张
设有两种系统聚类法A和B，他们在第i步的距
离矩阵分别为Ai和Bi（I=1，2，3…），若Ai>Bi ，
则称第一种方法A比第二种方法B使空间扩张，或
第二种方法比第一种方法浓缩。
3、方法的比较
D（短） 
D（长） 
当 0   1
当 0
D（平），D（重）  D（平）；
D（平）；
，D（变平） D（平）；
，D（变平） D（平）。
六、主要的步骤
1、选择变量
（1）和聚类分析的目的密切相关
（2）反映要分类变量的特征
（3）在不同研究对象上的值有明显的差异
（4）变量之间不能高度相关
2、计算相似性
相似性是聚类分析中的基本概念，他反映了研究
对象之间的亲疏程度，聚类分析就是根据对象之
间的相似性来分类的。有很多刻画相似性的测度
3、聚类
选定了聚类的变量，计算出样品或指标之间的相
似程度后，构成了一个相似程度的矩阵。这时主
要涉及两个问题：
（1）选择聚类的方法
（2）确定形成的类数
4、聚类结果的解释和证实

对聚类结果进行解释是希望对各个类的特征
进行准确的描述，给每类起一个合适的名称。这
一步可以借助各种描述性统计量进行分析，通常
的做法是计算各类在各聚类变量上的均值，对均
值进行比较，还可以解释各类产别的原因。

如果是变量聚类分析，聚类分析做完之后，各
类中仍有较多的指标。也就是说聚类分析并没有达
到降维的目的。这就需要在每类中选出一个代表指
标，具体做法是：假设某类中有 k 个指标，首先分
别 计 算 类 内 指 标 之 间 的 相 关 指
数  ij 2 (i  j ) (i  1,2,, k )，然后计算某个指标与类内
其他指标之间相关指数的平均数，即
Ri 
2

 ij
i j
k 1
取 xi 最大的 Ri 2 ，做为该类的代表。
例 某公司下属30个企业，公司为了考核下属企业的经
济效益，设计了8个指标。为了避免重复，需要对这8个指标
进行筛选，建立一个恰当的经济效益指标体系。通过计算30
个企业8个指标的相关系数距离，数据是1-r2。得如下表:
x1
x1
x2
x2
x3
x4
x5
x6
x7
x8
0
0.60
0
x3
x4
0.43
0.46
0
0.47
0.45
0.12
0
x5
x6
x7
x8
0.57
0.45
0.23
0.22
0
0.38
0.40
0.21
0.29
0.22
0
0.31
0.79
0.65
0.70
0.80
0.66
0
0.45
0.45
0.27
0.23
0.14
0.19
0.77
试用将它们聚类。
0
根据美国等20个国家和地区的信息基础设施
的发展状况进行分类。
Call—每千人拥有的电话线数；
move l—每千人户居民拥有的蜂窝移动电话数；
fee—高峰时期每三分钟国际电话的成本；
comp—每千人拥有的计算机数；
mips—每千人计算机功率（每秒百万指令）；
net—每千人互联网络户主数。
国家
meiguo
riben
deguo
ruidian
call
631.6
498.4
557.6
684.1
movel
161.9
143.2
70.60
281.8
fee
0.36
3.57
2.18
1.4
comp
403
176
199
246
mips
26073
10223
11571
16660
net
35.34
6.26
9.84
29.39
ruishi
644
498.4
93.5
147.5
1.98
2.5
234
284
13621
13578
22.68
13.49
xiongyali
469.4
434.5
81.9
138.6
92.2
174.9
169
262.2
56.1
73
16.3
8.20
9.8
5
6.5
49.4
3.68
3.36
3.02
1.4
2.61
5.12
3.68
2.66
119
99
19
31
31
24
40
68
6911
5795
876
1411
1751
1101
1796
3067
1.72
1.66
0.52
1.28
0.35
0.48
1.45
3.09
malaixiya
195.5
88.4
4.19
53
2734
1.25
taiguo
yindu
faguo
yingguo
78.6
13.6
559.1
521.10
27.8
0.30
42.9
122.5
4.95
6.28
1.27
0.98
22
2
201
248
1662
101
11702
14461
0.11
0.01
4.76
11.91
xinjiapo
taiwan
hanguo
baxi
zhili
moxige
eluosi
bolan
Pseudo F Statistic
120
110
100
90
80
70
60
50
40
30
20
10
0
1
2
3
4
5
6
7
8
9
10
11
Number of Clusters
12
13
14
15
16
17
18
19
§2 动态聚类
一、思想
系统聚类法是一种比较成功的聚类方法。然而当
样本点数量十分庞大时，则是一件非常繁重的工作，
且聚类的计算速度也比较慢。比如在市场抽样调查
中，有4万人就其对衣着的偏好作了回答，希望能迅
速将他们分为几类。这时，采用系统聚类法就很困
难，而动态聚类法就会显得方便，适用。
动态聚类解决的问题是：假如有个样本点，要把
它们分为类，使得每一类内的元素都是聚合的，并
且类与类之间还能很好地区别开。动态聚类使用于
大型数据。
Yes
选
择
凝
聚
点
分
类
是
否
合
理
分
类
修
改
分
类
No
分
类
结
束
用一个简单的例子来说明动态聚类法的工作
过程。例如我们要把图中的点分成两类。快速
聚类的步骤：
1、随机选取两个点 x (1) 和 x2(1)作为聚核。
1
2、对于任何点 x ，分别计算d ( x , x )和d ( x , x )
(1)
k
k
(1)
1
k
2
3、若d ( x , x )  d ( x , x )，则将 x 划为第一类，
(1)
k
1
(1)
k
2
k
否则划给第二类。于是得图（b）的两个类。
4、分别计算两个类的重心，则得 x 和 x ，以其为
(2)
(2)
1
2
新的聚核，对空间中的点进行重新分类，得到新分类。
（a）空间的群点
(c) 第一次分类
(b) 任取两个聚核
(d) 求各类中心
(e) 第二次分类
二、选择凝聚点和确定初始分类
凝聚点就是一批有代表性的点，是欲形成类的中
心。凝聚点的 选择直接决定初始分类，对分类结果也
有很大的影响，由于凝聚点 的不同选择，其最终分类
结果也将出现不同。故选择时要慎重．通 常选择凝聚
点的方法有：
(1) 人为选择，当人们对所欲分类的问题有一定了
解时，根据经验，预先确定分类个数和初始分类，并
从每一类中选择一个有代表性的样品作为凝聚点。
(2) 将数据人为地分为A类，计算每一类的重心，
就将这些重心作为凝聚点。
(3) 用密度法选择凝聚点：以某个正数d为半径，
以每个样品为球心，落在这个球内的样品数(不包括作为
球心的样品)就叫做这个样品的密度。计算所有样品点的
密度后，首先选择密度最大的样品作为第一凝聚点，并
且人为地确定一个正数D(一般D＞ d，常取D＝2d)。然
后选出次大密度的样品点，若它与第一个凝 聚点的距离
大于D，则将其作为第二个凝聚点；否则舍去这点，再
选密度次于它的样品。这样，按密度大小依次考查，直
至全部样品考查完毕为止．此方法中，d要给的合适，
太大了使凝聚点个数太 少，太小了使凝聚点个数太多。
(4) 人为地选择一正数d，首先以所有样品的均值
作为第一凝聚点。然后依次考察每个样品，若某样品
与已选定的凝聚点的距 离均大于d，该样品作为新的凝
聚点，否则考察下一个样品。
(5) 随机地选择，如果对样品的性质毫无所知，可
采用随机数表来选择，打算分几类就选几个凝聚点。
或者就用前A个样品作为凝聚点(假设分A类)。这方法
一般不提倡使用。
确定初始分类常用的方法有：
(1) 人为分类，凭经验将样品进行初步分类．
(2) 选择凝聚点后，每个样品按与其距离最近的凝聚点归类。
(3) 选择一批凝聚点后，每个凝聚点自成一类，将样品依次
归入其距离最近的凝聚点所在的类，并重新计算该类的重心，以
代替原来的凝聚点，再考虑下一个样品的归类，直至所有样品都
归类为止．
(4) 用某种聚类方法得到一个分类，这个分类就可作为初始
分类．当样本量大时，有时只用部分样品按某种聚类方法进行分
类．用每类重心作为凝聚点，再用(2)或(3)的方法对全部样品归类
后得到初始分类．
三、衡量聚类结果的合理性指标和算法终止的标准
定义
设 P 表示在第n次聚类后得到的第i类集
n
i
合，i  1,2,3,, k ，Ai(n ) 为第n次聚类所得到的聚核。
定义
k
un  

i 1xPi n
d 2 ( x, Ai( n ) )
1
 xl
Ai 
ni x P
j
j
l
j  1,2,
i
第j＋1步的新聚核。
k
若分类不合理时，un  
i 1xPi
d ( x, A ) 会很大，随着分
2

n
(n)
i
类的过程，逐渐下降，并趋于稳定。
定义
第i类中所有元素与其重心的距离的平方和:
D( Ain , Pi n )
k
un   
i 1xl Pi
xl Pi n
d 2 ( xl , Ain )
d ( xl , A ) 
2
n

(n)
i
k
n
 D( Ai
i 1
n
, Pi )
为所有K个类中所有元素与其重心的距离的平方和。
| un1  un |

算法终止的标准是
un1
 是事前给定的一个充分小量 。
五、动态聚类步骤为：
第一，选择若干个观测值点为“凝聚点”；
第二，可选择地，通过分配每个“凝聚点”最近的类里
来形成临时分类。每一次对一个观测值点进行归类，“凝
聚点”更新为这一类目前的均值；
第三，可选择地，通过分配每个“凝聚点”最近的类里
来形成临时分类。所有的观测值点分配完后，这些类的“
凝聚点”用临时类的均值代替。该步骤可以一直进行直到
“凝聚点”的改变很小或为零时止；
第四，最终的分类有分配每一个观测到最近的“凝聚点
”而形成。
例 我国经济发展的总目标是到2000年人民生活达到
小康标准，因此，了解各地区目前对小康生活质量
的实现程度。对各地区实现小康生活质量的状况进
行综合评价，对各级政府部门具有重要意义。数据
是1990年全国30个省在经济（jj）、教育（jy）、
健康（jk）和居住环境（jz）四个方面对小康标准
已经实现的程度，1表示已经达到或超过小康水平，
0表示低于或多或少刚达到温饱水平。希望利用该数
据对15个地区进行分类研究。
jj
jy
jk
jz
类别
距离
beijngs
0.7258
0.9413
1.0000
0.5000
1
0.29550
shangh
0.5346
0.9848
1.0000
0.5000
1
0.14909
tianjin
0.3246
0.9733
1.0000
0.5000
1
0.16173
henna
0.2301
0.4621
1.0000
1.0000
2
0.22252
zhej
0.5025
0.2374
1.0000
0.8882
2
0.34448
jilin
0.3446
0.7755
0.8280
0.5000
1
0.18212
helongj
0.2891
0.7835
0.8080
0.5000
1
0.22322
fujian
0.1406
0.3524
1.0000
0.7102
2
0.27468
guangx
0.0939
0.6498
0.4435
1.0000
2
0.51560
anhui
0.1104
0.0802
1.0000
0.9545
2
0.34050
ningxia
0.2708
0.3127
0.5425
0.9053
2
0.29445
hunan
0.0618
0.5687
0.4385
0.5000
3
0.41704
jiangxi
0.0549
0.3042
0.3520
0.6155
3
0.15540
Qinghai
0.0751
0.0118
0.0000
0.8258
3
0.37720
Guiz
0.0286
0.0600
0.0590
0.5000
3
0.25968
四、有序样本聚类法
（一）功能范畴与数据类型
有序样本聚类法又称为最优分段法。该方法是由费歇
在1958年提出的。它主要适用于样本由一个变量描述的情
况。或者将多变量综合成为一个变量来分析。
设
   1 , 2 ,, n  是样本点构成的集合，样本点 i在
函数 V ( )上的取值为  i 。若 V ( i )  V ( j ) ，则将视为一类。
不妨假设 v1   2     m 。要将 v1 , 2 ,, m 分为 K 类；
即P  ( P1 , P2 ,, Pk )
，分类时不能打乱样本点的顺序，即每
一类必须呈的  1 , 2 ,, j  形式，即有序样本聚类。
系统聚类开始n个样品各自自成一类，然后
逐步并类，直至所有的样品被聚为一类为止。而
有序聚类则相反，开始所有的样品为一类，然后
分为二类、三类等，直到分成n类。每次分类都
要求产生的离差平方和的增量最小。
例    1 , 2 ,, 4 
 1  V ( 1 )  1.1
 2  V ( 2 )  V ( 3 )  1.3
 3  V ( 4 )  3.1
这里n=4，m=3。若将其分为两类，其结果应该是
P  (v1 , v2 ), v3
对应中的点是( w1 , w2 , w2 ), w4 。
有序样本聚类法常常被用于系统的评估问题，被用来对
样本点进行分类划级。例如，十二个地区的经济发展指数，
排列出来以后，需要划分他们的等级。一种方法是按照行
政命令。规定三个经济发达地区，四个中等发达的地区，
三个一般地区，两个发展较差地区。
这种行政上的规定往往是不客观、不合理的。合理的分类
应该把发展情况最近似的地区划入同一类。这就是有序样
本聚类的工作思路。
（二）有序聚类的步骤
设有序样品x(1),x(2),…,x(n)。
1、定义类的直径
设某类G中包含的样品有
x
(i)
, x(i 1) ,, x(j) 
( j  i)
该类的均值向量为
j
1
 x(t)
XG 
j  i  1 t i

用D(i,j)表示这一类的直径，常用的直径有：
D(i, j )   (x(t)  XG )(x(t)  XG )
j
t i
 当是单变量的时，也可以定义直径为：
~
D(i, j )   | x(t)  X G |
j
t i
~
X G是中位数
2、定义分类的损失函数
用b(n,k)表示将n个有序的样品分为k类的
某种分法：
G1  i1 , i1  1,, i2  1
G1  i2 , i2  1,, i3  1


Gk  ik , ik  1,, n
定义这种分类法的损失函数为
Lb(n, k )   D(it , it 1  1)
k
t 1
当n和k固定时， L[b(n,k)]越小表示各类的离差平方
和越小，分类是合理的。因此要寻找一种分法b(n,k)，使
分类损失函数L[b(n,k)]达到最小。记该分法为p[n,k]。
3、L[b(n,k)]的递推公式

D(1, j  1)  D( j , n)
 L p(n,2)  min
2 j  n

D( j  1, k  1)  D( j, n)
 L p(n, k )  min
k  j n
以上的两个公式的含义是，如果要找到n个样品分为k
个类的最优分割，应建立在将j-1（j＝2,3,…,n)个样品分为
k-1类的最优分割的基础上。
4、最优解的求法
若分类数k是已知的，求分类法b(n,k)，使它
在损失函数意义下达到最小，其求法如下：
首先，找出分点jk，使
L[ P(n, k )]  L[ P( jk  1, k  1)]  D( jk , n)
于是得第k类 Gk   jk , jk  1,, n
然后，找出jk－1，使它满足
L[ P( jk 1 , k  1)]  L[ P( jk 1  1, k  2)]  D( jk 1 , jk  1)
于是得第k-1类 Gk 1   jk 1 , jk 1  1,, jk  1
类推。一次可以得到所有类G1，G2，…Gk，
这就是所求得最优解。
分析儿童的生长期。有如下的资料是1-11
岁的男孩平均每年的增重：
年龄
1
2
3
4
5
6
7
8
9
10
11
增加重量
（公斤）
9.3
1.8
1.9
1.7
1.5
1.3
1.4
2.0
1.9
2.3
2.1
问男孩的发育可分为几个阶段。
直径D(i,j)
1
1
0.444
2
0.598
3
0.709
4
0.819
5
0.943
6
1.078
7
1.225
8
1.412
9
1.638
10 1.915
2
3
4
5
6
7
8
9
10
0.005
0.02
0.051
0.104
0.179
0.287
0.427
0.613
0.839
0.005
0.021
0.053
0.104
0.181
0.287
0.433
0.614
0.005
0.021
0.051
0.105
0.181
0.292
0.434
0.005
0.02
0.055
0.106
0.187
0.294
0.005
0.024
0.056
0.111
0.187
0.008
0.023 0.005
0.057 0.021 0.006
0.107 0.051 0.021 0.005
样品数
3
4
5
6
7
8
9
10
11
2类
0.005(2)
0.020(2)
0.051(2)
0.104(2)
0.179(2)
0.287(2)
0.427(2)
0.613(2)
0.839(2)
3类
4类
最小误差函数表
5类 6类 7类
0.005(4)
0.010(4)
0.026(5)
0.041(5)
0.075(5)
0.107(6)
0.160(7)
0.211(7)
0.005(5)
0.010(5)
0.015(6)
0.033(7)
0.046(8)
0.062(8)
0.092(8)
0.005(6)
0.010(6) 0.005(7)
0.015(8) 0.010(8) 0.005(8)
0.021(8) 0.015(8) 0.010(8) 0.005(9)
0.037(8) 0.021(10) 0.015(10) 0.010(10) 0.005(10)
0.051(10) 0.026(10) 0.021(11) 0.015(11) 0.010(11) 0.005(11)
8类
9类
10类
损失函数L[p(n,k)]随k变化趋势图
1
损失函数
0.8
0.6
0.4
0.2
0
1
2
3
4
5
类数
6
7
8
9
分类数
2
3
4
5
6
7
8
9
10
误差函数
0.8392
0.21069
0.09199
0.05102
0.02586
0.02055
0.01523
0.01016
0.00508
最优分割结果
1,2-11
1,2-6,7-11
1,2-4,5-7,8-11
1,2-4,5-7,8-9,10-11
1,2-3,4-5,6-7,8-9,10-11
1,2-3,4-5,6-7,8-9,10,11
1,2-3,4,5,6-7,8-9,10,11
1,2-3,4,5,6,7,8-9,10,11
1,2-3,4,5,6,7,8,9,10,11
其他聚类方法
两维图论聚类
由于生态区划的划分既要求各分区内部的
自然、经济条件具有相对一致性，保持空间
上的连通性和行政区界的完整性，如果采用
普通的聚类方法，就很难使分类结果满足分
区的原则，因此，需引人图论的空间表示法
及相关算法。
0-1变量聚类
以0-1变量为原始变量进行聚类分析。
模糊聚类
在自然科学或社会科学研究中，存在着许多定义不很严
格或者说具有模糊性的概念。这里所谓的模糊性，主要是指
客观事物的差异在中间过渡中的不分明性，如某一生态条件
对某种害虫、某种作物的存活或适应性可以评价为“有利、
比较有利、不那么有利、不利”；灾害性霜冻气候对农业产
量的影响程度为“较重、严重、很严重”，等等。这些通常
是本来就属于模糊的概念，为处理分析这些“模糊”概念的
数据，便产生了模糊集合论。在侧重于应用的模糊数学分析
中，经常应用到聚类分析、模式识别和综合评判等方法。
二阶段聚类：（资料来自spss.in)
1.目的：在样本容量太大时替代系统聚类。（具有系
统聚类和K－均值聚类的优点）
2.特征：
1）变量类型：连续型，离散型（或数量型，质量型）
2）通过对观测的“两层”聚类，有效的降低了样本容
量太大对层次聚类的限制；同时，进一步提高了聚类的效
果，使组间的差异较组内差异更加显著。
3.样品与样品间，类与样品间，类与类间的相互联系程
度的量度：
1）对数似然距离——实质是概率，注意：这里有
关变量的原始假定: （1）连续型变量服从正态分布
（2）类别型变量服从多项分布
（3）变量相互独立
2）欧氏距离——这里程序默认：类与类，类与样
品间的距离采用重心距离。
4.具体程序步骤：
1）将样本观测归入数量较多的子类。方法：有条
件的经修改的层次聚类——
(1)
程序对样本进行异端值（或异端类）的
识别，以避免异端值对聚类树的影响（可使聚类树
“体积”过大——类内的差异不合理的扩大，对聚类
效果造成影响）
(2)
层次聚类
进行经修改（详情请阅读“BIRCH”）的
(3) 由给定阈值（或初始类数）的条件下，完成初始聚类。
注意：一般来说，子类越多，二阶段聚类越准确，但同时对计
算机硬件的要求越高，故应合理确定初始类的数量。(这里程
序建议使用赤池或schwarts准则）
2)在给定阈值（pass）（也可根据其他标准确定类数）的
条件下，对初始类（如果在（1）中采用了对异端值、类的处
理，这里的初始类将不包括它们）进行层次聚类——类间距离
将采用重心距离（程序默认）。
3）确定合理的类数.
愿同学们喜欢数学建模
预祝在数学建模竞赛中
取得好成绩