数学建模与计算

数学建模与计算

1 多元线性回归

1.1 回归系数的计算

假设变量y 与p 个变量x 1, x 2, … , x p 之间存在以下线性关系：

y = b 0 + b 1x 1 + b 2x 2 + … + b p x p (+ ε) (1.1)

这些变量的n 组观测值如表1.1所示。

将这些数值代入(1.1)得

y 1 = b 0 + b 1x 11 + b 2x 12 + … + b p x 1p y 2 = b 0 + b 1x 21 + b 2x 22 + … + b p x 2p ……

y n = b 0 + b 1x n 1 + b 2x n 2 + … + b p x np (1.2)

这是一个关于b 0, b 1, …,b p （待估参数）的线性方程组。一般来说方程个数n 远大于变量个数p , 是一个不相容的线性方程组。但是可以求出b 0, b 1, …,b p 的一组数值代入(1.2)式右边，使得等号两边的数值尽可能接近。最常用的便是最小二乘法，即

Min Q (b 0, b 1, …,b p ) =

∑=----n

i ip p i i x b x b b y 1

2110)...(

令

Y = ????

??

? ??n y y y 21, X =

??

?

?

?

??

??np n n p p x x x x x x x x x 212222*********, β =

??????

?

??p b b b 10, 得

Min Q (β) = (Y - X β)T (Y - X β) = Y T Y - 2Y T X β + βT X T X β

这是一个关于β的凸二次函数。令

β

β??)

(Q = -2 X T Y + 2 X T X β = 0， 得到

X T X β = X T Y .

此式称为正规方程组或法方程组。如果X T X 可逆，则

β

?= (X T X )-1X T Y . (1.3) 其中β?= (0?b ,1?b ,…, p

b ?)T 。 定义1. n 元实值函数f (x ) = f (x 1, x 2, … , x n )的n 个偏导数构成的向量称为梯度函数，记为?f (x )或

x

x f ??)

(, 即 x x f ??)( = T

n

x x f x x f x x f ))(,,)(,)((21?????? . 结论1.1. 设c 和x 是n 维数向量，那么x

x c T ??)

(= c .

结论1.2. 设A 是n ? n 对称矩阵，x 是n 维向量，那么x

Ax x T ??)

(= 2Ax 。

结论1.3. X T X 可逆的必要充分条件是X 的列线性无关。

1.2 回归方程的检验

定义

总平方和SS tot = ∑-=n

i i y y 1

2

)(,

回归平方和SS reg = ∑-=n i i y y

12

)?(, 残差平方和SS err = ∑-=n

i i i y

y 1

2

)?(, 其中y =∑=n i i y 1/ n , i y

? = 0?b +1?b x i 1 + … +p b ?x ip 。另外 R 2

=

tot

reg SS SS ,

调节的R 2 = 1 - (1 - R 2)

1

1

---p n n .

通常，R 2或调节的R 2大于0.95时可认为回归方程成立。

1.3 可以用线性回归方法求解的非线性模型

(1) 二次函数

y = c 0 + c 1x + c 2x 2.

将x 2看作一个新的变量，这是一个具有2个自变量的线性模型。

设x 的观测值为x 1, x 2, … , x n , y 的观测值为y 1, y 2, … , y n , 则公式(1.3)中的

X = ??????

?

?

?222

2211111n n x x

x x x x , Y = ????

??

?

??n y y y 21, β = ????

? ??210c c c

(2) C-D 函数

Y = AL αK β

其中L 是劳动力的投入，K 是资金投入，Y 是产出，A 称为技术系数。A , α, β是待估参数。两边取对数得

ln Y = ln A + αln L + βln K

(3) 施工区交通事故模型[2]

ln f = α0 + α1ln L + α2D

1 + α3Q 1+ α4U 其中f 是施工区交通事故频率，L 是施工区长度，D 是施工期，Q 是交通流量, U 是道路类型。

(4) 单隧道地表沉降槽的高斯曲线模型

S = S max 2

ax e

其中S 是x 处的沉降深度，S max 和a 是待估参数。两边取对数得

ln S = ln S max + ax 2.

1.4 不可以用线性回归方法求解的非线性模型

(1) 逻辑增长曲线预测模型

bt

t ae k

y -+=

1

其中y t 是第t 期某经济指标的数值，k , a , b 是待估参数。

(2) 双隧道地表沉降槽的高斯曲线模型[3]

S = s 12

1x a e

+ s 22

2)(u x a e

-

其中S 是x 处的沉降深度，u 是两隧道中心的距离，s 1, a 1, s 2, a 2是待估参数。

以上两个模型的参数可用非线性最小二乘法求取。假设变量y 自变量x (∈R p )之间存在以下关系:

y = f (x ; β)(+ ε)

其中β ∈ R m 是待估参数。

给定因变量与自变量的若干组观测值(y i , x i )，非线性最小二乘的数学模型如下：

∑-=i

i i x f y Q 2)),((min ββ

这是一个无约束极值问题，由于形式特殊，可用Levenberg-Marquardt(LM)方法求解，见[4]p262-266.

习题

(1) 证明结论1.1, 1.2和1.3。

(2) 编一个计算程序计算(X T X)-1。

参考文献

[1] 盛骤，谢式千，潘承毅. 概率论与数理统计，高等教育出版社，1997.

[2] Qiang Meng, Jinxian Weng, Xiaobo Qu. A probabilistic quantitative risk assessment model for the long-term work zone crashes, Accident Analysis and Prevention, 2010, (42) 1866-1877.

[3] 马灵. 基于数据挖掘的隧道施工地表沉降规律研究，华中科技大学博士学位论文，2013.

[4] Jorge Nocedal, Stephen J Wright. Numerical Optimization数值最优化，科学出版社，2006.

2 神经网络

2.1 BP神经网络

回归分析是寻求因变量y与自变量x (∈R p)之间函数关系的一种统计分析方法。这种方法需要给定两者之间函数关系的形式y = f(x; β), 利用因变量与自变量的若干组观测值(y i, x i)确定参数β(∈R m)，一般用最小二乘法。y = f(x; β)是x的线性函数时称为线性最小二乘问题，是x的非线性函数时称为非线性最小二乘问题。

对于许多实际问题，因变量与自变量之间的关系很复杂，难以确定它们之间函数关系的形式。（人工）神经网络不要求明显地给出因变量与自变量之间函数关系的表达式。神经网络是由一些单元（神经元）按一定方式连接而成的网络，其形式多种多样，最常用的是分层前向神经网络，其第一层是输入层，最后一层是输出层，中间一般有1或2层，称为隐层。输入层和输出层单元的个数由需要解决的问题本身决定，隐层单元的个数由实验确定，使得网络的输出尽可能与目标值接近。图2.1是一个具有三个输入两个输出一个隐层的神经网络，其中的圆圈表示单元，与圆圈相连的箭线表示其输入或输出。

图2.1 神经网络的例子

神经网络的每个单元对应一个数字，称为输出。输入层单元的输入相当于自变量的观测值。隐层和输出层每个单元的输入是前一层各单元输出的加权和，称为净输入。将净输入代入一个称为激励函数或传递函数的一元函数，得到的数值便是该单元的输出。为了防止各单元的净输入不至于过大或过小，有时需加上一个数值，称为调整量。

神经网络的主要工作是不断调整权值，称为训练，使得其输出尽可能接近理想输出（相当于回归分析中因变量的观测值）。调整权值的方法很多，最有名的是分层神经网络中使用的误差反向传播（Back Propagation, BP ）算法。

在分层神经网络中，隐层和输出层每个单元的输出递归地确定。设神经网络有p 个输入，q 个输出，分L 层。激励函数为g (?)。用w ij l 表示第l - 1层单元i 指向第l 层单元j 的连接权。在给定p 个输入y 11, y 21, … , y p 1和所有连接权的前提下，其他单元的输出由左到右由上到下依次计算。

第l 层单元j 的净输入

h j l =

∑-i

l ij l i w y 1 第l 层单元j 的输出

y j l = g (h j l ), l = 2, … , L .

对隐层和输出层每个单元由左到右进行的以上运算称为前向传播。

以图2.1网络部分单元为例说明前向传播计算方法。隐层第1个单元的净输入

h 12 = y 11w 112 + y 21w 212 + y 31w 312.

该单元的输出

y 12 = g(h 12).

输出层第1个单元的净输入

h 13 = y 12 w 113+ y 22w 213 + y 32w 313 + y 42w 413.

该单元的输出

y 13 = g(h 13).

计算涉及的数据见图2.2。

2

图2.2 前向传播（部分）示意图

一般来说，图2.1隐层各单元的净输入和输出为

∑==3

1

212i ij i j

w y h , )(2

2j j h g y =, j = 1, 2, 3, 4. (2.1)

图2.1输出层各单元的净输入和输出为

∑==4

1

323i ij i j

w y h ，)(3

3j j h g y =，j = 1, 2. (2.2) 激励函数有许多，例如（[1]p124）

(1) 线性函数

g (x ) = cx .

(2) 符号函数

g (x ) = sgn(x ) = ?

??<-≥.0,1,

0,1x x

(3) Sigmoid 函数

g (x ) =

x

e -+11.

(4) 双曲函数

g (x ) = x

x

e e --+-11.

神经网络每层用到的激励函数可能不一样（[2]p132），根据实验效果确定。

设神经网络最后一层单元j 的目标值为d j L ，为了使得各单元的输出y j L 尽可能接近d j L ，极小化误差平方和：

∑==-=q j L j L j w

y d E 12)(21min ∑=-q

j L j L j h g d 12

))((21 (2.3)

以获取连接权w ij l 的数值。E 的表达式称为成本函数或业绩函数，w 是所有连接权构成的向

量。由于

h j L =∑-k

L

kj L k w y 1,

按照复合函数求导法，E 对第L - 1层单元i 与第L 层单元j 的连接权w ij L 的偏导数

L

ij

w E

??= j h E ??L ij L j w h ??= -(d j L - y j L ) g ’(h j L ) y i

L -1. (2.4) 其中

L j

h E ??=-(d j L - y j L ) g ’(h j L )。在(2.4)中令

δj L = g ’(h j L ) (d j L - y j L )，

得到

L

ij

w E

??= -y i L -1δj L .

于是，按梯度法确定的w ij L 的改变量为

?w ij L = η y i L -1δj L ,

其中η是0与1之间的数字，称为学习率。

现考虑第l - 1层单元i 与第l 层单元j 的连接权w ij l 的改变量, l = L - 1, … , 2。利用复合函数求导法求成本函数E 对w ij l 的偏导数，按梯度法w ij L 的改变量为

?w ij l = η y i l -1δj l , l = l - 1, (2)

其中

δj l = g ’(h j l )∑++i

l i l ji w 1

1δ.

以上运算是由右到左进行的，称为反向传播。

以上调整连接权用到的方法就是求解无约束极值的梯度法。梯度法仅利用目标函数在当前解处关于各个变量的一阶偏导数确定调整量。所以一阶偏导数的计算是关键。下面以图2.1网络为例详细说明计算方法。其成本函数为

E (w ) =

])()[(2

12323223131y d y d -+-. 此式隐含有3?4 + 4?2 = 20个变量w ij l 。

首先求E 对w ij 3的偏导数。由于E 是y j 3的函数，y j 3是h j 3的函数, h j 3是w ij 3的函数,

3ij

w E

??= -(d j 3 - y j 3)g ’(h j 3)y i 2. 令

δj 3 = (d j 3 - y j 3)g ’(h j 3),

得

3ij

w E

??= -δj 3y i 2. 于是

?w ij 3 = ηδj 3y i 2,

w ij 3 := w ij 3 + ?w ij 3, i = 1, 2, 3, 4, j = 1, 2.

再求E 对w ij 2的偏导数。由于E 是y 13和y 23的函数，y 13和y 23分别是h 13和h 23的函数, h 13

和h 23是y j 2的函数, y j 2是h j 2的函数, h j 2是w ij 2的函数,

2ij

w E

??= -(d 13 - y 13)g ’(h 13) w j 13g ’(h j 2)y i 1 - (d 23 - y 23)g ’(h 23) w j 23 g ’(h j 2) y i 1. (注：将(2.2)前面的式子改为∑

==4

1

31231j j j w y h ,∑==4

1

3

223

2

j j j w y h ) 由于

(d j 3 - y j 3)g ’(h j 3) = δj 3, j = 1, 2,

得

2ij

w E ??= -δ13 w j 13g ’(h j 2)y i 1 - δ23 w j 23 g ’(h j 2) y i 1

= -(δ13 w j 13 + δ23 w j 23)g ’(h j 2)y i 1.

令

δj 2 = (δ13 w j 13 + δ23 w j 23)g ’(h j 2),

得

2ij

w E

??= -δj 2 y i 1. 于是

?w ij 2 = ηδj 2y i 1,

w ij 2 := w ij 2 + ?w ij 2, i = 1, 2, 3, j = 1, 2, 3, 4.

上面的δj l 是第l 层第j 个单元的“误差”，l = 2, 3。图2.3表示计算δ12要用到w 113, δ13, w 123, δ23。

2

图2.3 误差反向传播（部分）示意图

现将BP 算法的计算步骤归纳如下[3]。

(1) 确定分层神经网络的结构，随机地令每个连接权为很小数值。 (2) 选择输入信号y i 1, i = 1, 2, … , p 。 (3) 沿网络前向传播信号。 (4) 计算输出层各单元的误差：

δj L = g ’(h j L )(d j L - y j L ), j = 1, 2 ,… , q .

(5) 由右到左计算隐层各单元的误差：

δj l = g ’(h j l )∑++i

l i l ji w 1

1δ, l = L - 1, … , 2.

(6) 更新连接权

w ij l → w ij l + ?w ij l , l = 2, … , L 其中?w ij l =η y i l -1δj l ，η是很小正数。

(7) 返回(2)，直到成本函数值E 足够小。

2.2 径向基函数网络

径向基函数(radial basis function, RBF)网络一般为三层：输入层，以非线性RBF 函数为激励函数的隐含层，以及线性输出层。

RBF 函数是形如Φ(x ; c ) = Φ(|| x - c ||)的函数，其中x , c ∈ R n , || ? ||是欧氏范数，c 称为中心。在一元函数的情形，该函数图形关于x = c 对称。

设输入层有n 个单元，隐含层有m 个单元，网络仅一个输出。则输出

y = f (x ) =

∑=-m

i i i c x Φw 1

||)(||.

其中x 是n 个输入构成的向量。

RBF 函数网络是要通过不断调整隐含层各单元RBF 函数的中心c i 以及其与输出单元的连接权w i ，使得网络的输出尽可能接近目标值。中心的调整和连接权的调整是分开进行的。计算步骤大致如下：

步骤1 初始步骤

从训练样本集中随机抽取m 个样本作为初始中心。 步骤2 调整中心 随机抽取一个样本，计算隐含层各单元的输出并调整中心。此步骤反复进行直到隐含层产生L (> m )组输出。

步骤3 调整连接权

利用最小二乘法计算连接权，然后返回步骤2，直到达到一定精度为止。 在RBF 函数网络中，最常用的RBF 函数是高斯函数

Φ(|| x - c ||) = }2||||exp{2

2

σc x --

. 给定整数K (> m )和足够大的整数L 。RBF 函数网络的计算流程如下(参考[2]p159)： (1) 从训练样本集中随机抽取m 个样本，记为c i (i = 1, 2, … ,m )，作为初始中心。令l = 1, k = 1。

(2) 从训练样本集中随机抽取一个样本x , 计算隐含层m 个单元的输出

?ki = }2||||exp{2

i

i c x σ--, i = 1, 2, … ,m , 其中σi =

m

S 2max

, S max 是m 个中心的最大距离。 设与样本x 对应的输出的目标值为d k 。 (3) 调整中心

从m 个中心中找出距离x 最近的中心：

i k = arg i

m in || x - c i ||,

并令

c i = ???≠==-+.

,,...,2,1,;),(k i k i i i i m i c i i c x c η

其中η ∈ (0, 1)。

k := k + 1。若k ≤ K , 转(2)。否则 (4) 调整连接权 由最小二乘模型

∑=-K

k k k w

d w 1

2)(min ?

求出新的连接权w = (w 1, w 2, … ,w m )T ，其中?k = (?k 1, ?k 2,…, ?km ), k = 1, 2, … ,K 。

(5) l := l + 1。若l ≤ L , 令k = 1, 转(2), 否则停止。

参考文献

[1] 吴今培，孙得山. 现代数据分析，机械工业出版社，2006 [2] 张得丰. MATLAB 神经网络编程，化学工业出版社，2011

[3] Anil K. Jain, Jianchang Mao. Artificial Neural Networks: A Tutorial, IEEE, 1996

3 主成分分析

3.1 主成分分析法原理[1]

主成分分析法利用p 个可观测的变量X 1, X 2,…, X p 构造m 个变量Y 1, Y 2,…, Y m (m ≤ p )，使得Y i (i = 1, 2,…, m )是X 1, X 2,…, X p 的线性组合，其组合系数构成的向量是单位向量(长度等于1)，其中Y 1的方差最大，Y 2的方差第2大且Y 2与Y 1不相关，Y 3的方差第3大且Y 3与Y 2和Y 1都不相关，以此类推。这样得到的Y 1, Y 2,…, Y m 称为X 1, X 2,…, X p 的第1个，第2个, … , 第m 个主成分。

考虑X 1, X 2,…, X p 的线性组合

Y = w 1X 1 + w 2X 2 + … + w p X p

其中w 12 + w 22 + … + w p 2 = 1。设E (X i ) = x i , 则Y 的方差为

D (Y ) =

E [Y - E (Y )]2

= E [(w 1X 1 + w 2X 2 + … + w p X p ) - (w 1x 1 + w 2x 2 + … + w p x p )]2 = E [(X 1 - x 1)w 1 + (X 2 - x 2)w 2 + … + (X p - x p )w p ]2

=

∑∑==--p i p

i j i j j i i w w x X x X E 11)])([( =

∑∑==p i p

i j i ij w w 11

σ

= w T Gw

其中，G = (σij )p ?p , 称为协方差矩阵，σij = COV (X i , X j ) = E [(X i - x i )(X j - x j )]。

现在要求一个单位向量w = (w 1, w 2,…, w p )使得D (Y ) = w T Gw 达到最大，即求解

Max w T Gw

s.t. w T w = 1. (3.1) 为此构造Lagrange 函数

L (w , λ) = w T Gw - λ (w T w - 1).

令

w

w L ??)

,(λ= 2Gw - 2λw = 0 得

Gw = λw .

可见要求的组合系数向量w 是协方差矩阵G 的特征向量。在上式两边左乘以w T 得到

w T Gw = λw T w .

当w T w = 1时λ = w T Gw 。所以模型(3.1)的目标函数值就是协方差矩阵G 的特征值。

用λ1表示G 的最大特征值，w 1 = (w 11, w 12, … , w 1p )T 表示相应的单位特征向量，则第1个主成分

Y 1 = w 11X 1 + w 12X 2 + … +w 1p X p .

第2个主成分是以下模型的解([2]p348)：

Max w T Gw

s.t. w T w = 1

w 1T w = 0. (3.2)

设其单位特征向量是w 2 = (w 21, w 22, … , w 2p ), 那么第2个主成分为

Y 2 = w 21X 1 + w 22 X 2 + … +w 2p X p .

值得注意的是协方差矩阵G 的特征值和相应特征向量通常不是用(3.1)和(3.2)等优化模型求得的，而是用QR 等方法求出来的。

3.2 主成分分析法计算步骤

设X 1, X 2,…, X p 的n 组观测值为

表3.1 p 个变量的观测值

步骤1. 计算样本协方差矩阵G 计算

∑==n

j ji i x n x 1

1, i = 1, 2,…, p

和

D = ??

??

?

?

?

?

?---------p np n n p p p p x x x x x x x x x x x x x x x x x x

2211222211212

12111 则G =

D D n T 1

1

-. 步骤2. 计算G 的特征值和相应特征向量

设计算得到的p 个特征值λ1 ≥ λ2 ≥ … ≥ λp , 相应单位特征向量为

w i = (w i 1, w i 2, … , w ip ), i = 1, 2, … ,p .

步骤3. 确定主成分 计算方差贡献率

∑

=p

j j

i

1λλ(i =1, 2, … , p )和累计方差贡献率

∑

=p j j

11

λλ+

∑

=p

j j

12

λλ+ … +

∑

=p

j j

k

1λλ.

当前k 个方差贡献率的和达到80%以上时，就得到k 个主成分

Y i = w i1X1 + w i2 X2 + … + w ip X p, i = 1, 2, … , k.

习题

1) A是n阶方阵，其特征值和特征向量是如何定义的？当A是实对称矩阵时，其特征值和特征向量有何特点？

2) X和Y是随机变量，它们的协方差和相关系数是如何定义的?

3) 设X的观测值为x1, x2, … , x n, Y的观测值为y1, y2, … , y n。写出样本协方差和样本相关系数的计算公式。

4) p个随机变量X1, X2, … , X p的协方差矩阵和相关系数矩阵是如何定义的?

参考文献

[1] 刘顺忠. 管理统计学和SAS软件应用，武汉大学出版社，2006.

[2] 蒋尔雄，高坤敏，吴景琨. 线性代数，人民教育出版社，1979.

4几种启发式优化算法

最优化分为无约束最优化和约束最优化两类，无约束最优化模型为

Min f(x) = f(x1, x2, … , x n)

约束最优化模型为

Min f(x) = f(x1, x2, … , x n)

s.t. x∈X?R n

其中x = (x1, x2, … , x n)T, f(x)称为目标函数, X称为可行域，可行域中的点称为可行解。Max f(x)与Min -f(x)等价。

常用的启发式(heuristic)算法有遗传算法、粒子群优化和模拟退火等等。

4.1遗传算法

遗传算法(Genetic Algorithm, GA)是John Holland等人于1970年代中期提出的一种模仿生物进化过程计算方法。

遗传算法的基本术语有染色体(chromosome)，基因(gene)，种群(population)，种群大小(population size)，适应度函数(fitness function)，交叉(crossover, recombination)，变异(mutation)和选择(selection)等。

染色体是可行解的编码，一般用由0和1组成的字符串表示。字符串中的每一位叫基因。种群是若干染色体构成的集合，种群中的染色体叫个体(individual)，种群大小则是种群中个体的数目。适应度函数就是目标函数，但在遗传算法中一般要求极大化。交叉和变异是对字符串进行操作的两种方式，用来改变种群中某些个体。交叉是对种群中两个个体进行的操作，以产生两个新的个体。前两个叫父母(parents)，后两个叫后代(offspring)。变异是对个体的某个位进行的取反操作。选择是从种群中随机挑选若干个体，再用以前找到的最好染色体取代其中适应值最低的个体，组成新种群的过程，每个个体被选中的可能性与它的适应值成正比。

例如，求解以下问题(GAlgorithm.pdf)：

Max f(x) = x2

s.t. x∈ {0, 1, … , 31}.

由于有32 = 25个可行解，用5位的0-1字符串(染色体)表示这32个数值：

00000 ? 0, 00001 ? 1, …. , 11111 ? 31.

假如我们决定将种群大小定为4，并且随机地从32个染色体中选择出4个构成一个初

始种群：01101, 11000, 01000, 10011。

为了计算每个个体的适应度，先对它们进行解码：

01101 → 13, 11000 → 24, 01000 → 8, 10011 → 19.

再计算f (x ) = x 2：

13 → 169, 24 → 576, 8→ 64, 19 → 361.

交叉和变异也是随机进行的。交叉分为一点交叉，两点交叉，均匀交叉等方式。设父母为01101和11000。一点交叉将它们在同样的位置截成2段，例如01101和11000，交叉结果为01100和11001。两点交叉将它们截成3段，交换中间一段。均匀交叉则将两个字符串均匀地截为多段，间隔地交换各段。

变异是对个体的某个位进行的取反操作，例如01101 → 00101。

对初始种群和经过交叉和变异形成的新种群，还要对其进行选择(selection, recombination)操作。其中一种为赌轮选择法([1]p167)。设种群大小为n , 这n 个个体的适应值为f 1, f 2, … , f n (如果有负数，每个加上同一个正数使它们成为正数)。计算步骤如下：

(1) 计算每个个体被选取的概率

∑

==

n

k k

i

i f f P 1, i = 1, 2, … , n

和累计概率

∑==

i

j j i P q 1

, i = 1, 2, … , n .

(2) 产生区间[0, 1]随机数r ，若r < q 1，则第1个个体入选，否则选择第i 个个体 (2 ≤ i

≤ n )，它满足q i -1 < r < q i 。

(3) 重复(2)直到有n 个个体为止。 得到n 个新的个体后，计算每个个体的适应值。用迄今为止适应值最高的个体替代其中适应值最低的个体。如果当前种群适应值最高的个体比迄今为止适应值最高的个体还高，则以其作为迄今为止适应值最高的个体。

遗传算法的计算步骤如下： (1) 如果目标是求极小，将其改为求极大。如果适应度有负数，则加上一个适当的正数，使其为正。

(2) 将问题的可行解转化为染色体(字符串)，设字符串长度为L 。 (3) 给定种群大小n ，交叉概率p c 和变异概率p m 。给定迭代次数N 。 (4) 随机生成n 个个体构成初始种群。令k = 0。

(5) 计算n 个个体的适应值f i 和总适应值F =

∑=n

i i f 1。

(6) 选择操作，计算每个个体的选择概率F f P i i /=和累计概率∑==i

j j i P q 1，并以赌轮法选择n 个个体。再用以前找到的最好染色体取代其中适应值最低的个体。 (7) 交叉操作

1) 依次对每个个体产生[0, 1]间随机数r ，若r < p c ，则该个体参入交叉操作，选出2个为止；

2) 对每一对个体，随机生成2到L 间一个整数，在该位置进行交叉操作。

(8) 变异操作

对每个个体中的每一位产生[0, 1]间的随机数r ，若r < p m 则在该位取反。 (9) 令k = k + 1, 若k > N 则停止，否则转(5)。

用遗传算法求函数极值时自变量的数值只能是有限个离散值，而且其个数需等于2n , n 是二进制数的长度。当某变量x 在区间[a , b ]内变化时可将此区间分为2n - 1等份，用2n 个分区间的端点作为x 在区间[a , b ]的代表。如果要求计算误差不超过ε，n 应满足

1

2--n a

b ≤ ε, 由此得

n ≥

2

ln )

1ln(

+-ε

a

b .

其中 δ =

1

2--n a

b 是分区间的长度。 例如a = -1, b = 2, ε = 10-6,

δ =2

ln )1103ln(6+?=21.5165.

可取n = 22。

4.2 粒子群优化

粒子群优化(Particle Swarm Optimization, PSO)是Kennedy 和Eberhart 于1990年代中期提出的一种模拟鸟群飞行行为的一种计算方法。

这种方法使用粒子(Particle)、位置(position)、速度(velocity)、粒子群(swarm)等概念。粒子和位置是可行解，速度是粒子移动的方向和大小，粒子群是若干个粒子的集合。

粒子群优化的计算步骤如下[2]： (1) 初始步骤

随机产生n 个粒子和速度：

x 0i = x min + rand (x max - x min ), v 0i =

t

x x rand x ?-+)

(min max min , i = 1, 2, … ,n .

其中, x 0i 和v 0i 分别是第i 个粒子的初始位置和速度，x min 和x max 分别是粒子各分量的最小值和最大值构成的向量，rand 是0和1之间的一个随机数，?t 是一个正数，一般就令?t = 1。

(2) 更新每个粒子的速度

v k +1i = wv k i + c 1rand

t x p i k i ?-+ c 2rand t

x p i

k

g k ?-, i = 1, 2, … ,n . 其中, w 称为惯性因子(inertia factor)，在0.4和1.4之间；c 1称为自信度(self confidence)，在1.5和2之间，c 2称为群信度(swarm confidence)，在2和2.5之间；p i 是到目前为止第i 个粒子的最好位置(适应值最小的位置)；p k g 是当前群中适应值最小的粒子。

(3) 更新每个粒子的位置

x k +1i = x k i + v k +1i ?t , i = 1, 2, … ,n .

(4) 返回步骤(2)。

以上过程反复进行，直到连续几次迭代(如10次迭代)粒子群的最好适应值变化不超过ε (如10-5)时为止。

文献[2]就8例子进行试验，比较遗传算法和粒子群优化的优劣。这8例子有无约束极值问题也有约束极值问题，前2个问题是

求Rosenbroch 函数的极小值：

Min f (x 1, x 2) = 100(x 2 - x 12)2 + (1 - x 1)2,

每个变量的下限为-5，上限为5。此问题的最优解为(x 1, x 2) = (1, 1)。

求Eggerate 函数的极小值：

Min f (x 1, x 2) = x 12 + x 22 + 25(sin 2x 1 + sin 2x 2)2,

每个变量的下限为-2π，上限为2π。此问题的最优解为(x 1, x 2) = (0, 0)。

其余6个是工程技术方面的问题。

试验中，种群和粒子群的大小n = 40。在粒子群优化中，w = 0.5, c 1 = 1.5, c 2 = 2.5。试验结果是粒子群优化的效果要好一些。

4.3 模拟退火

模拟退火(Simulated Annealing, SA) 由N. Metropolis 等人于1953年提出，是模拟物理过程的一种计算方法。

设问题为求极小值。模拟退火的计算步骤如下([3]p15)：

令T min 等于一个适当正数，给定正数T < T min , 并随机产生一个初始解x 。 (1) 产生x 的一个临近点y 。

(2) 如果f (y ) < f (x )，令x := y 。否则 (3) 计算

P y = ))

()(exp(T

x f y f --

并随机产生在0与1之间的一个数r 。如果P y > r ，令x := y 。

(4) 稍微减少T 值，如果T < T min ，返回步骤(1), 否则停止。

注：步骤(2)执行的是下降，但步骤(3)可能上升。由以上公式可以看出，如果f (y )比f (x )大很多或者T 值较小，计算得到的P y 值小，上升的可能性较小。参数T 类似温度。温度高时固体内部粒子呈现出一种无序状态，温度逐渐降低时，状态趋于稳定。

参考文献

[1] 吴今培，孙德山. 现代数据分析，机械工业出版社，2006.

[2] Rania Hassan, Babak Cohanim, Olivier de Week. A Comparison of Particle Swarm Optimization and the Genetic Algorithm （网上下载）

[3] Wesley E. Snyder, Hairong Qi. 机器视觉教程，机械工业初步社，2005

5 粗糙集

5.1 粗糙集的一些基本概念

一个知识表达系统或信息系统S 可以表示为有序四元组

S = {U , R , V , f }

式中，U = {x 1, x 2, … , x n }是由n 个(研究)对象构成的论域, R = C ? D 为属性集，其中C 是条件属性集，D 是决策属性集，V 是由映射f : U ? R → V 确定的属性值的集合。系统S 的表格

形式称为决策表，例如表5-1反映的医疗信息是一个决策表。

对于属性集R 的某个子集R 1，一般R 1 ? C ，如果对象x i 和x j 的这些属性值都相同，就说x i 和x j 是在R 1下是不可分辨的，并将它们归为一类。这样便可得到论域U 的一种分类(划分)，记为U / R 1, 其中的每个子集叫R 1基本集，记为1][R x 。给定X ? U ，所有包含于X 的基本集的并称为下近似集，记为R _(X ), 即

R _(X ) = {x ∈ U | 1][R x ? X }.

所有与X 的交不为空集的基本集的并称为上近似集，记为)(X R -, 即

)(X R -= {x ∈ U | 1][R x ? X ≠ ?}.

正域Pos (X ) = R _(X ), 负域Neg (X ) = U - R _(X ), 边界域Bnd (X ) = )(X R -- R _(X )。

如果上近似集等于下近似集，即)(X R -= R _(X ), 称X 是可定义的，否则称X 是不可定义的。不可定义集有以下4种类型[6]：

(1) 如果R _(X ) ≠ ?且)(X R -≠ U , X 称为是粗糙定义的； (2) 如果R _(X ) ≠ ?且)(X R -= U , X 称为是外部不可定义的； (3) 如果R _(X ) = ?且)(X R -≠ U , X 称为是内部不可定义的； (4) 如果R _(X ) ≠ ?且)(X R -= U , X 称为是完全不可定义的。 集合X 的粗糙程度用粗糙度)(1X R α表示，

)(1X R α =

))

(())_((X R Card X R Card -

其中Card (?)表示一个集合的基数，即集合中所含元素的数目。显然0 ≤)(1X R α≤ 1。如果

)(1X R α = 1，集合X 相对于R 1是可定义的，如果)(1X R α< 1，集合X 相对于R 1是粗糙的。 例5.1 考虑表5.1。对于属性子集R 1 = {头疼, 肌肉疼}, 计算样本子集X = {x 1, x 2, x 5}的

下近似集、上近似集、正域、边界域和粗糙度。

解 按属性子集R 1划分的论域为

U / R 1 = {{x 1, x 2, x 3}, {x 4, x 6}, {x 5}}.

样本子集X 与3个基本集的关系是

{x 1, x 2, x 5} ? {x 1, x 2, x 3} = {x 1, x 2} ≠ ?， {x 1, x 2, x 5} ? {x 4, x 6} = ?， {x 1, x 2, x 5} ? {x 5} = {x 5} ≠ ?.

由于3个基本集中仅{x 5} ? X , R _(X ) = {x 5}。由于

{x 1, x 2, x 5} ? {x 1, x 2, x 3} ≠ ?，{x 1, x 2, x 5} ? {x 5} ≠ ?，

)(X R -= {x 1, x 2, x 3}?{x 5} = {x 1, x 2, x 3, x 5}。另外正域

Pos (X ) = R _(X ) = {x 5},

边界域

Bnd (X ) = )(X R -- R _(X ) = {x 1, x 2, x 3}.

粗糙度

)(1X R α =

))

(())_((X R Card X R Card -=

4

1

= 0.25. 例5.2 表5-2是6个考生情况调查表。计算这些考生按升学情况分类的分类精度。

解 按高考成绩将考生分类:

U / R = {{3, 5}, {2}, {1, 6}, {4}} = {Y 1, Y 2, Y 3, Y 4},

按升学情况将考生分类:

U / X = {{2, 3, 5, 6}, {1, 4}} = {X 1, X 2}.

(1) X = X 1 = {2, 3, 5, 6} (升学考生集合)时，因为Y 1 ? X 1, Y 2 ? X 1，下近似集

R _(X 1) = Y 1 ? Y 2 = {2, 3, 5}.

因为Y 1 ? X 1 ≠ ?, Y 2 ? X 1 ≠ ?, Y 3 ? X 1 ≠ ?, 上近似集

)(1X R -= Y 1 ? Y 2 ? Y 3 = {1, 2, 3, 5, 6}.

于是边界域

Bnd (X 1) = )(1X R -- R _(X 1) = {1, 6},

粗糙度

)(1X R α= Card (R _(X 1)) / ))((1X R Card - = 3 / 5.

即考上大学的分类精度为3/5。

(2) X = X 2 = {1, 4} (未升学考生集合)时，因为Y 4 ? X 2, 下近似集

R _(X 2) = Y 4 = {4}.

因为Y 3 ? X 2 ≠ ?, Y 4 ? X 2 ≠ ?, 上近似集

)(2X R -= Y 3 ? Y 4 = {1, 4, 6}.

于是边界域

Bnd (X 2) = )(2X R -- R _(X 2) = {1, 6},

粗糙度

)(2X R α= Card (R _(X 2)) / ))((2X R Card - = 1 / 3.

即未考上大学的分类精度为１/ 3。

5.2 决策表的约简

决策表的约简是指去掉一些冗余属性，使得按剩余属性对论域U 所做的划分与原属性集R 对U 所做的划分不变。设R i , R j , R k 是属性集R 的子集，如果按{R i , R j }和按{R i , R k }对论域U 作的分类与按R 对U 进行的分类相同，{R i , R j }和{R i , R k }称为R 的简约。此处{R i , R j }?{R i , R k } = R i ≠ ?，称R i 为核，并记为Core (R ) = R i 。

例5.3 设有知识库K = (U , R ), 其中论域U = {x 1, x 2, … , x 8}, 属性集R = {R 1, R 2, R 3}。假设论域U 按属性子集R 1, R 2, R 3的划分分别为

U / R 1 = {{x 1, x 4, x 5, x 6}, {x 2, x 3}, {x 7, x 8}}, U / R 2 = {{x 1, x 2, x 5}, {x 4, x 6, x 7}, {x 3, x 8}}, U / R 3 = {{x 1, x 2, x 5}, {x 4, x 6}, {x 3, x 7}, {x 8}}.

求R 的简约和核。

解 由于x 1和x 5都在3种划分的某个基本集中，x 4和x 6也都在3种划分的某个基本集中，可见按R 划分的论域为

U / R = {{x 1, x 5}, {x 2}, {x 3}, {x 4, x 6}, {x 7}, {x 8}}.

另外，

U / (R - R 1) = U / {R 2, R 3} = {{x 1, x 2, x 5}, {x 3}, {x 4, x 6}, {x 7}, {x 8}} ≠ U / R , U / (R - R 2) = U / {R 1, R 3} = {{x 1, x 5}, {x 2}, {x 3}, {x 4, x 6}, {x 7}, {x 8}} = U / R , U / (R - R 3) = U / {R 1, R 2} = {{x 1, x 5}, {x 2}, {x 3}, {x 4, x 6}, {x 7}, {x 8}} = U / R .

因此，{R 1, R 3}和{R 1, R 2}都是R 的简约。核Core (R ) = {R 1, R 3} ? {R 1, R 2} = R 1。

利用分辨矩阵和分辨函数对决策表进行约简十分方便。设论域U = {x 1, x 2, … , x n }，条件属性集C = {c 1, c 2,…, c m }，决策属性集D = {d }。分辨矩阵是一个n ? n 对称矩阵，其i 行第j 列的元素m ij 是x i 和x j 的具有不同属性值的属性的集合(析取)。分辨函数是分辨矩阵中所有非空元素的合取，其简化表示可确定决策表的简约及核。

例5.3 设有信息系统S = (U , R ), U = {x 1, x 2, … , x 6}, R = {a , b , c , d }。有关数据如表5.3所示。利用分辨矩阵和分辨函数求简约及核。

解 依次比较表5.3的两行，得到的分辨矩阵见表5.4。

分辨函数f S 是以上所有非空元素(集合)的合取。利用布尔公式x ∧ (x ∨ y ) = x 可将其化为

f S = b ∧ (a ∨ d )

再利用公式x ∧ (y ∨ z ) = (x ∧ y ) ∨ (x ∧ z )得

f S = (a ∧ b ) ∨ (b ∧ d ) = ab ∨ bd .

因此该信息系统有{a , b }和{b , d }两个简约，核是{b }。两个约简数据表分别为表5.5和表5.6。

表5.5 简约表1 表5.6 简约表2

5.3 决策规则

设某问题的约简决策表有m 个条件属性u 1, u 2, … , u m ，一个决策属性v 。如果(u 1, u 2, … , u m ) = (1u ,2u , … , m u )时，恒有v =v , 就说(1u ,2u , … , m u ) →v 是一个决策规则。 例5.4 某问题的约简决策表见表5.7（[1]p21），其中b 和c 是条件属性，e 是决策属性。求决策规则。

表5.7 某问题的约简决策表

解 由此表可看出以下一些决策规则（[1]p22）：

(b , 0)∧(c , 2) → (e , 1), (b , 2)∧(c , 0) → (e , 2), (b , 2)∧(c , 2) → (e , 0), (b , 1)∧(c , 0) → (e , 2), (b , 1)∧(c , 1) → (e , 2), (b , 1)∧(c , 2) → (e , 1).

若干术语

粗糙集(Rough Set)，对象(object)，论域(universe)，条件属性(condition attribute)，决策属性(decision attribute), 分类(classification)，划分(partition)，基数(cardinality)，不可分辨的(indiscernible), 可定义的(definable)，不可定义的(un(definable), 约简(reduce)，简约(reduct)。

参考文献

[1] 吴今培，孙德山. 现代数据分析，机械工业出版社，2006.

6 支持向量机

6.1 可被平面分隔的两类点

假设R d 空间有两个点集，可以用一个平面分隔开来。求一个平面将这两类点分开，并且距离这两类点尽可能远。

设有n 个点x i , i = 1, 2, … ,n , 如果x i 属于第一个点集，记y i = +1, 如果x i 属于第二个点集，记y i = -1. 问题的模型如下:

,min

ββ2

1βT β s.t. y i (x i T β + β0) ≥ 1, i = 1, 2, … , n . (6.1)

其中，β ∈ R d , β0 ∈ R 。

以上模型是一个凸二次规划。设最优解为β*, β0*，满足x i T β* + β0* = ±1的x i 称为支持向量。

对于某个x ∈ R d , 如果x T β* + β0* ≥ 1, 则x 属于第一类点; 如果-(x T β* + β0*) ≥ 1, 则x 属于第二类点, 否则不能判断x 的归属。

6.2 不可被平面完全分隔的两类点