概率统计模型(1)

概率统计模型

自然界中存在两种现象：确定性现象和不确定性现象．同一实验或者试验在不同次重复中，可能出现不同的结果的现象称为随机现象．随机现象的结果尽管是不确定的，但是，同一随机现象的多次重复却表现出某种规律性，即同一事件在不同次试验或者实验中出现的概率是确定的、唯一的．因此，随机现象中包含确定性现象．对随机现象的研究可以通过对随机现象的某些事件的发生概率来研究．变量之间也存在两种关系：确定性关系和不确定性关系．确定性关系：可用一个表达式确切描述，如圆的面积与半径之间的关系．描述确定性关系的数学模型有函数，微分方程，差分方程等．不确定性关系：不可用一个表达式确切描述，如人的体重与身高等．不确定性关系在现实生活中大量存在，即使许多看来是确定性关系的变量，在实际中也会受到各种不同随机因素的影响而变得不确定，确定性关系只是它们的一种近似，如自然科学的很多规律．

本章主要介绍利用概率统计知识分析随机现象和随机数据，建立随机模型，求解随机模型，并对得到的结果进行分析，最后运用于实际．第一节介绍几个直接利用概率知识的建模问题，如赌博问题，巴拿赫(Banach)火柴盒问题，信与信封的配对问题，切割机的收益问题；第二节回归分析模型，主要介绍施肥效果分析问题；第三节判别分析模型，主要介绍螨虫分类问题；第四节时间序列分析，主要介绍Chesapeake 海湾的收成预测问题；第五节随机模拟模型，主要介绍利用随机模拟方法产生随机数据及模拟随机现象的方法；第六节排队论模型，主要介绍用排队论的方法分析，处理等候问题．

通过以上这些模型和方法的学习，使读者了解和掌握一些处理随机问题的一般思想和方法，如果读者想进一步学习和了解随机数学的专业理论与方法，可阅读随机数学的一些分支的专门著作，如：随机过程，时间序列分析，回归分析，多元统计分析等．

§4-1 几个直接利用概率知识的建模问题

对随机现象的研究可以通过对随机现象的某些事件的发生概率来研究．本

节就来介绍几个概率模型，主要利用的基本知识就是古典概率模型的概率计算及其相关问题，随机变量的概率分布及其计算．可以参看任意一本大学理工科的《概率论与数理统计》教科书[7]，也可以参考周义仓、赫孝良两位老师编写的教科书[6]． 问题描述

问题1：赌博问题

均匀正方体骰子的六个面分别编号1，2，3，4，5，6．现将一对骰子抛掷6次以决定胜负，请问将赌注押在“出现两个１点”和“完全不出现两个１点”哪个更有利？

问题2：巴拿赫(Banach)火柴盒问题

波兰数学家巴拿赫随身带着两盒火柴，分别放在两个衣袋里，每盒有n 根火柴．使用时，每次随机地从其中一盒中取出一根．试求他将一盒火柴用完时，另一盒剩余火柴根数的分布律．

问题3：信与信封的配对问题

某人给它的N 个朋友写信，写好后，分别将这些信装入N 个信封中，并在信封上随机、不重复地写上N 个收信人的地址．问他一个都没写正确和恰有r 个写正确的概率各是多少？

问题4：切割机的收益问题[3]

一台线切割机把金属线切割成规定的长度．由于切割机的某种不准确性，切割线的长度X 可以看作是在区间[11.5,12.5]上的均匀分布的随机变量．规定的长度是12cm ．如果11.7≤X

元的损失丢弃．试计算：如果切割N 段金属线，那么，请估计平均每根金属线为老板贡献的利润是多少？ 问题求解

1. 问题１的求解

问题１是一个古典概率模型的概率计算问题．解决这样的问题的关键就是事件的表示．为此，我们令A i k 分别表示第i 次抛掷骰子时第k 枚骰子（k =1, 2）

出现１点的事件．那么，在第i 次抛掷中，两枚骰子都出现１点的事件A i 表示为

A i =A i 1A i 2 (4.1.1)

而6次抛掷中至少出现一次两个１点的事件B 可以表示为

B = (A i 1A i 2) (4.1.2)

i =16

这样，事件B 的对立事件是

=∏(+1i

i =16

2i

)=∏ (4.1.3)

i i =1

6

所以

P (B )=1-P () (4.1.4)

由于事件A i -1, A i -2, i =1,2,3,4,5,6相互独立，于是有

P (B )=1-P ()=1-∏P (i ) (4.1.5)

i =16

而

55⎛5⎫35

P (i )=P ()+P ()-P ()P ()=+- ⎪=

66⎝6⎭36

i

2i

1i

2i

2

所以

⎛35⎫

P (B )=1-P ()=1-∏P (i )=1- ⎪=0.1555 (4.1.6)

⎝36⎭i =1

6

6

这样，出现两个一点的概率是0.1555，大大小于完全不出现两个一点的概率0.8445．因此，应将赌注押在＂完全不出现两个一点＂上．

2. 问题2的求解

设巴拿赫总共取出的火柴根数为Z ，而分别来自于两个火柴盒，设从左右口袋的两个火柴盒中分别取出的火柴根数分别是X , Y ，于是

Z =X +Y (4.1.7)

而用U 表示总共剩余的火柴根数，于是

(U =k )=(X =k , Y =n ) (Y =k , X =n ) (4.1.8) 那么，我们要计算的随机变量U 的分布列．设巴拿赫发现左口袋火柴刚好取完时，右口袋里还剩Y =k 根火柴，因此，右口袋已经被取了n -k 根．这样，当巴拿赫首次发现左口袋没有火柴时，已经进行了n -k +n =2n -k 次随机试验．在这2n -k 次试验中，事件A ，即火柴取自左口袋出现了n 次，事件，即火柴取自右口袋出现了n -k 次．对右口袋先取完，我们又类似的讨论．这样，这个

问题实际上是一个二项分布的概率计算问题．于是

P (U =k )

=P (X =k , Y =n )+P (Y =k , X =n )

=P (X =k |Y =n )P (Y =n )+P (Y =k |X =n )P (X =n ) (4.1.9) =C =C

n 2n -k

⎛1⎫ 1-⎪⎝2⎭⎛1⎫ ⎪⎝2⎭

n -k

⎛1⎫1n ⎪+C 2n -k ⎝2⎭2

n

⎛1⎫ ⎪⎝2⎭

n -k

⎛1⎫1 1-⎪⎝2⎭2

n

2n -k

n 2n -k

, k =0,1, 2, , n

3. 问题3的求解

经过分析，问题相当于将N 封写好的信放到写着正确地址的信封．问题要求，计算所有的信都没有正确放到该放的信封的事件的概率，以及计算恰有r 封信正确放到该放的信封的概率．

这是一个古典概型问题．我们分别用, B r 表示没有正确放到该放的信封的事件以及恰有r 封信正确放到该放的信封的事件．用A i 表示第i 封信能正确放对信封的事件，i =1, 2, , N ，那么

=12 N (4.1.10)

这里特别要注意：1, , N 不是相互独立的，而事件A i 1A i 2 A i r i r +1 i N 是互不

相容（i 1, i 2, , i N 是1, 2, , N 的一个排列）的．因此，不能利用下式计算概率P ()

P ()=P (1)P (2) P (N ) (4.1.11)

⎛N ⎫

P ()=1-P (A )=1-P A i ⎪ (4.1.12)

⎝i =1⎭

但是，注意到公式

而

N N

⎛N ⎫N N -1

P A i ⎪=∑P (A i )-∑P (A i A j )+∑P (A i A j A k )+ +(-1)P (A 1 A N )

i , j =1, i

N N

1

=N ⋅-∑P (A i )P (A j |A i )+∑P (A i )P (A j |A i )P (A k |A i A j )+

N i , j =1, i

+(-1)

2

=1-C N

N

N -1

P (A 1) P (A N |A 1A 2... A N -1)

111111N -1131⋅+C N ⋅⋅+ +(-1)⋅ N N -1N N -1N -2N N -11

k -1

=∑

k =1

(-1)

k !

(4.1.13)

所以，由(4.1.12)，有

P ()=∑

k =0N

(-1)

k !

k

(4.1.1４)

用C r 表示恰好指定的r 封信装对信封，则由乘法原理，B r 中的样本点数为（这里n (B r ) 和n (C r ) 分别是事件B r , C r 的基本事件个数，或称样本点数）

r

n (B r )=C N n (C r ) (4.1.15) 而

P (B r )=

n (B r )N !

r

=C N

n (C r )N !

k

(4.1.16)

根据前面的分析和结论，有

P (C r )=∑

k =0N -r

(-1)

k !

k

(4.1.17)

而由古典概率的计算公式，有

P (C r )=∑

k =0N -r

(-1)

k !

=

n (C r )

N -r ! (-1)

k !

k

(4.1.18)

于是，得到

n (C r )=(N -r )! ∑

k =0k

N -r

(4.1.19)

1N -r (-1)P (B r )=∑, k =0,1, 2, , N (4.1.20) r ! k =0k !

4. 问题4的求解

我们只要知道了在三类区间的线段的数目，就可以计算出总的收益．设长度在区间11.7≤X

N =N p +N g +N l (4.1.21)

如果总利润是I ，那么平均每根金属线的利润为

N p N g N I

w ==0.25+0.10-0.02l (4.1.22)

N N N N

N N N

p , g , l 分别是随机变量X 落在如上三个区间的频率，而频率具

N N N

有稳定性，当N 充分大时，频率近似等于相应的概率值，即

N p N N g N N l N

1

dx =0.5

11.712.5-11.512.51

≈P (X ≥12.2)=⎰dx =0.3

12.212.5-11.511.71

≈P (X

11.512.5-11.5≈P (11.7≤X

12.2

所以，平均来说，单根金属线的利润为

w =0.25⨯0.5+0.10⨯0.3-0.02⨯0.2=0.151（元） (4.1.23)

研究性问题

4-2-1 供电问题：设某车间有200台车床相互独立地工作，由于经常需要维修、测量、调换刀具、变换位置等种种原因要停车．若每台车床有60%的时间在开动，而每台车床在开动时要耗电1KW ，问应供给这个车间多少电力才能保证在８h 生产中大约仅有0.5min 因电力不足而影响生产？

4-2-2 钓鱼问题：为了估计湖中鱼的数量，先从湖中钓出r 条鱼做上记号，并放回湖中过一段时间后再从湖中钓出S 条鱼，结果发现其中有x 条鱼标有记号．问应该如何估计湖中鱼的数量N ．

§4-2 农作物施肥量与产量的关系

问题描述

某地区农作物生长所需的营养素主要是氮（N ）、磷（P ）、钾（K ），农作物研究所在该地区对土豆与生菜做了一定数量的实验，实验数据如表4.2.1，其中：ha 表示公顷，t 表示吨，kg 表示公斤．当一个营养素的施肥量变化时，将另二个营养素的施肥量保持在第七水平，如对土豆关于N 的施肥量做实验时，P 与K 的施肥量分别取196kg/ha(第七水平) 与372kg/ha(第七水平) ．

表4.2.1 施肥量与产量实验数据

土豆

N

施肥量

(kg/ha) 0 34

产量 (t/ha) 15.18 21.36

施肥量 (kg/ha) 0 24

P

产量 (t/ha) 33.46 32.47

施肥量 (kg/ha) 0 47

K

产量 (t/ha) 18.98 27.35

67 101 135 202 259 336 404 471

N

施肥量 (kg/ha) 0 28 56 84 112 168 224 280 336 392

25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75

49 73 98 147 196 245 294 342

P

36.06 37.96 41.04 40.09 41.26 42.17 40.36 42.73

96 140 186 279 372 465 558 651

K

34.86 38.52 38.44 37.73 38.43 43.87 42.77 46.22

产量 (t/ha) 11.02 12.70 14.56 16.27 17.25 22.59 21.63 19.34 16.12 14.11

施肥量 (kg/ha) 0 49 98 147 195 294 391 489 587 685

产量 (t/ha) 6.39 9.48 12.46 14.33 17.10 21.94 22.64 21.34 22.07 24.53

施肥量 (kg/ha) 0 47 93 140 185 279 372 465 558 651

产量 (t/ha) 15.75 16.76 16.89 16.24 17.56 19.20 17.97 15.84 20.11 19.40

试建立模型分析施肥量与产量的关系，并对所得结果从应用价值与如何改进等方面作出分析． 问题分析

农作物的产量与施肥量之间存在密切的关系，但很难用一个确定的函数关系来表达， 故可考虑用回归分析方法来研究其相关关系，建立回归方程近似描述产量与施肥量之间的相关关系. 模型假设

1. 实验中，只考虑施肥量对农作物产量的影响，其它因素：如温度，湿度，其它微量元素的含量，均处于相同水平，不预考虑.

2. 各次实验相互独立，结果互不影响，观测误差独立同分布，服从N (0, σ2), σ>0，N ，P ，K 的用量可精确控制，误差忽略不计.

变量及符号说明

n : 实验总次数，本问题中为10.

Q 1i : 对土豆而言，第i 次实验的产量，i =1, 2, , n Q 2i : 对生菜而言，第i 次实验的产量，i =1, 2, , n

Q 1Ni : 对土豆而言，与N 1i 对应的第i 次实验的产量，i =1, 2, , n Q 1Pi : 对土豆而言，与P 1i 对应的第i 次实验的产量，i =1, 2, , n

Q 1Ki : 对土豆而言，与K 1i 对应的第i 次实验的产量，i =1, 2, , n Q 2Ni : 对生菜而言，与N 2i 对应的第i 次实验的产量，i =1, 2, , n Q 2Pi : 对生菜而言，与P 2i 对应的第i 次实验的产量，i =1, 2, , n Q 2Ki : 对生菜而言，与K 2i 对应的第i 次实验的产量，i =1, 2, , n N 1i : 对土豆而言，第i 次实验的N 的用量，i =1, 2, , n N 2i : 对生菜而言，第i 次实验的N 的用量，i =1, 2, , n P 1i : 对土豆而言， 第i 次实验的P 的用量，i =1, 2, , n P 2i : 对生菜而言， 第i 次实验的P 的用量，i =1, 2, , n K 1i : 对土豆而言， 第i 次实验的K 的用量，i =1, 2, , n K 2i : 对生菜而言， 第i 次实验的K 的用量，i =1, 2, , n .

模型建立

1. 先对实验数据，作出散点图，直观分析产量与施肥量的变化趋势及关系. 从散点图来看，三种营养素的施肥量与产量之间存在非线性关系，尤其，氮肥的施用量与产量之间存在明显的二次关系，故可考虑建立三种营养素的施肥量与产量之间的一元二次回归模型．

2. 三种营养素的施肥量与产量之间的一元二次回归模型

Q 1Ni =a 10+a 11N 1i +a 12N 12i +ε1Ni , i =1,2, , n

2

Q 1Pi =b 10+b 11P i =1,2, , n 1i +b 12P 1i +ε1Pi ,

Q 1Ki =c 10+c 11K 1i +c 12K 12i +ε1Ki , i =1,2, n Q 2Ni =a 20+a 21N 1i +a 22N 12i +ε2Ni , i =1,2, , n

2Q 2Pi =b 20+b 21P i =1,2, , n 1i +b 22P 1i +ε2Pi ,

Q 2Ki =c 20+c 21K 1i +c 22K 12i +ε2Ki , i =1,2, , n

对上述模型，由已知实验数据，利用Mathematica 软件编程计算可得回归方程．但

是，考虑到作物的产量是各种营养素综合作用的结果，而以上建立的仅仅是一元回归模型，故须对模型进行改进．

3. 包含所有变量的全回归模型

2

Q 1i =a 0+a N N 1i +a P P 1i +a K K 1i +a NN N 1i +a NP N 1i P 1i

2

+a NK N 1i K 1i +a PP P 1i +a PK P 1i K 1i +a KK K 1i +ε1i

由全回归模型的求解结果(如表4.2.7) 及残差可看出，残差均匀分布在零点两侧，

无系统偏差，模型基本合适．但注意到，作物产量受各种营养素的影响不是同样的，且营养素两两之间的交互作用对产量的影响也不是同等的，故需对变量进行选择，进行逐步回归．

4. 逐步回归模型

利用MA TLAB 中的逐步回归函数stepwise 对变量进行逐步回归，回归结果表明：

① 对于土豆，首先进入模型的是N 与K 的交互作用项，其次是NN 项； ② 对于生菜，首先进入模型的是P ，其次是NN 项． 模型求解

对以上三个模型的求解，采用MA TLAB 软件进行．结果如下：

1. 一元回归模型的结果及分析

对土豆而言，N 的施肥量与产量的回归方程系数：

表4.2.2

常数项 14.7416

一次项 0.1972

二次项

-0.0003

对土豆而言，P 的施肥量与产量的回归方程系数：

表4.2.3

常数项 32.9161

一次项 0.0719

二次项

-0.00013783

对土豆而言，K 的施肥量与产量的回归方程系数：

表4.2.4

常数项 24.4144

一次项 0.0749752

二次项

-7*10^(-5)

对生菜而言，N 的施肥量与产量的回归方程系数：

表4.2.5

常数项 79.2501

一次项 3.516472

二次项

-0.0106883

对生菜而言，P 的施肥量与产量的回归方程系数：

表4.2.6

常数项 6.87795

一次项 0.0606347

二次项

-5.5*10^(-5)

对生菜而言，K 的施肥量与产量的回归方程系数：

表4.2.7

常数项 16.2329

一次项 0.00511548

二次项

-7.2*10^(-7)

以上一元回归模型结果表明：二次项系数较小且为负值，说明产量先随施

肥量增加而增加，达到一个峰值，然后，随施肥量增加而下降．说明，在一定范围内，施肥量对产量有促进作用，这对我们在生产管理中，科学、有效、经济地确定施肥量具有指导意义．

2. 对土豆的全回归模型的结果及分析

表4.2.8

一次项

常数项

N 0.0749752

P 0.0265478

K 0.0284431

交互作用项 NP 0.000222494

NK 0.000173897

0 P K

NN

二次项

PP

KK

15.2093

-000

325779

-0.0

0017 1209

-0.0

00067809

结果表明：一次项系数由大到小依次是N ，K ，P ，交互作用项依次是NP ，NK ，说明我们在生产管理中，不但要重视每中肥料的单独作用，还要充分重视肥料间的交互作用，这样才能在生产中充分发挥肥料对产量的促进作用．

3. 逐步回归模型的结果

对于土豆，首先进入模型的是N 与K 的交互作用项，其次是NN 项; 对于生菜，首先进入模型的是P ，其次是NN 项．回归结果表明，对土豆等块茎类作物，NK 的交互作用对作物的生长起显著作用，对生菜等叶类作物，P 的作用非常显著，其次，N 的作用对各种作物都是重要的．得到的结果符合作物栽培学原理与实际经验． 研究性问题

以上是从产量的角度考虑其与施肥量的关系．对此问题，还可以从经济学的角度考虑以下问题：

1. 研究产量与肥料用量的变化关系，确定各种肥料的边际用量；

2. 考虑到各种肥料的成本不同，为了达到最大效益，确定各种肥料用量的 最佳组合．

§4-3 AF 螨虫和APF 螨虫的区分问题

现有9只AF 螨虫和6只APF 螨虫的触角长与翼长数据： AF ：（1.24，1.72），（1.36，1.74），（1.38，1.64），（1.38，1.82），（1.38，1.90），

问题描述

（1.40，1.70），（1.48，1.82），（1.54，1.82），（1.56，2.08）．

APF ：（1.14，1.78），（1.18，1.96），（1.20，1.86），（1.26，2.00），（1.28，2.00），（1.30，1.96）．

对以上数据，制定一种方法正确区分螨虫；依据确立的方法，判别新样品（1.4，1.80），（1.28，1.84），（1.40，2.04）的归属；若AF 是宝贵的益虫，APF 是某疾病的载体，是否修改分类方法． 问题分析

此问题属于判别分析问题，即根据样本的指标（螨虫的触角长与翼长），建立判别规则，来判断样本来自哪个总体（AF ，APF ）．判别分析的一般模型可这样描述：设有k 个总体G 1, G 2, , G k ，它们的分布分别是F 1(x ), F 2(x ), , F k (x )，均为p 维分布，制定判别规则，对给定的新样品，确定它来自哪个总体．判别分析的方法有很多，如距离判别，Bayes 判别，Fisher 判别等．这里，我们采用距离判别． 模型假设

1. 两种螨虫的触角长与翼长服从二维正态分布N 2(μ1, ∑1)，N 2(μ2, ∑2)，μ1≠μ2, ∑1≠∑2；

2. 判别时仅考虑触角长与翼长两项指标，不考虑其它指标． 模型建立

设AF 螨虫为总体G 1，APF 螨虫为总体G 2，G 1 N 2(μ1, ∑1)，G 2 N 2(μ2, ∑2)．

1. 首先对两总体的均值进行显著性检验，即检验：μ1=μ2，当其有显著性差异时再进行判别．

2. 给出样品X 到总体G i 的距离（这里采用马氏距离）

d i 2=(X -μi )' ∑i -1(X -μi ), i =1,2

3. 建立判别函数及判别规则 判别函数为

2

W (X )=d 2-d 12

判别规则为

⎧⎪⎨⎪⎩

若W (X )>0, 则X ∈G 1

若W (X )

模型求解

μ1=μ2，1. 首先，对两总体的均值进行显著性检验，即检验：利用MATLAB

软件统计工具箱中的kstest2函数检验两总体分布是否相同，利用ttest2检验均

值是否相同，检验结果表明：两总体分布相同，均值存在显著性差异，故可继续进行判别．

2. 利用已知样本数据，计算判别函数值．

由于两总体均值与方差未知，采用极大似然估计，即

ˆ=ˆi =i , ∑μi

1(i )

L xx , i =1, 2 n i -1

最终的判别函数为：W (x 1, x 2)=-2.9358+29.1128x 1-190.293x 2．

对最初的两类样本，代入，回判结果如下：

表4.3.1

样本序号

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

样本值 （1.24，1.72） （1.36，1.74） （1.38，1.64） （1.38，1.82） （1.38，1.90） （1.40，1.70） （1.48，1.82） （1.54，1.82） （1.56，2.08） （1.14，1.78） （1.18，1.96） （1.20，1.86） （1.26，2.00） （1.28，2.00） （1.30，1.96）

原属类别 判别函数值 AF AF AF AF AF AF AF AF AF APF APF APF APF APF APF

0.433676 3.54663 6.03181 2.60654 1.08419 5.47231 5.51782 7.26459 2.89922 -3.61936

判定类别 AF AF AF AF AF AF AF AF AF APF APF APF APF APF APF

-5.88012

-3.39494 -4.31227 -3.73002 -2.38659

对新样品的判别结果：

表4.3.2

样本序号 16 17 18

样本值 （1.4，1.80） （1.28，1.84） （1.40，2.04）

判别函数值 3.56938

判定类别 AF APF APF

-0.685328 -0.997652

结果分析

对制定的判别函数及判别规则，用已知的经验样本进行计算，验证，结果表明，回判正确率100%，判别规则及方法有效． 研究性问题

在判别分析中，应考虑误判损失，若AF 是宝贵的益虫，APF 是某疾病的载体，则本属于APF 而误判为AF 的损失要大于本属于AF 而误判为APF 的损失，则应提高进入AF 的阀值（即判别样本落入某一类的判别函数临界值，如 以上判别样本属于AF 的阀值为0）．

§4-4 Chesapeake 海湾的收成预测问题

时间序列分析的方法来源人们对生产实践中所产生的历史数据的分析．人

们期望通过这些数据获得对未来某个较近时间的数据的估计．一般地，我们所得到的数据可以写为下面的数据序列

(x 1, y 1), (x 2, y 2), , (x n , y n ) （4.4.1）

这里，x i 是n 维向量，y i 实值标量．

我们可以这样想象：历史数据（4.1.1）是按照某种具有固定生产程序的机

f (t |x ) （4.4.2）

器所产生，对于同一个x ，所对应的y 是按照一个条件分布密度函数 产生的，因而y 的期望值为

因此，从理论上，我们要寻找的依赖关系应当是这个函数关系(4.4.3)．这个函数关系称为回归函数．

我们的目的是借助于概率统计的方法给出实值变量y 与n 维向量x 之间的回归函数或者估计，并且给出这个函数或者估计的误差限．寻找这样的函数关系或者估计的方法是较多的．时间序列预测[4,8]的一些方法：如回归估计、平稳

y =⎰tf (t |x )dt =ϕ(x ) （4.4.3）

时间序列的滑动平均、自回归、自回归滑动平均模型、Markov 链等可以用来处理预测预报问题，也可以使用现代基于支持向量机[5-6]的非参数统计的线性回归或非线性回归的方法．

本节将利用一般的基于最小二乘法的参数回归估计方法、以及基于支持向量机回归的非参数统计学习等方法来解决Chesapeake 海湾的收成预测问题，并主要介绍非平稳时间序列的预测问题．在这里我们不过多地拘泥于理论的陈述，具体的细节，请读者参阅相关资料． 问题描述

1992年《每日评论》（Daily Press）报告了过去50年中收集到的Chesapeake 海湾海产品收成方面的数据．我们将考察几种场合，并使用Chesapeake 海湾的商贸行业提供的如下数据：

（a ）收获蓝鱼的观测数据表4.4.1，（b ）收获蓝蟹的观测数据表4.4.1，回答下面两个问题：

表4.4.1 Chesapeake 海湾海产品收成方面的数据[11]

年 1940 1945 1950 1955 1960 1965

蓝鱼（磅） 15000 15000 250000 275000 270000 280000

蓝蟹（磅） 100000 850000 1330000 2500000 3000000 3700000

年 1970 1975 1980 1985 1990

蓝鱼（磅） 290000 650000 1200000 1500000 2750000

蓝蟹（磅） 4400000 4660000 4800000 4420000 5000000

问题1：请预测1995年收获的蓝鱼磅数； 问题2：请预测1995年收获的蓝蟹磅数．

（注 1磅＝453.6g ．）

问题分析

直观上，这不是一个平稳时间序列．因此，我们不能采用处理平稳时间序列的模型[8]进行预测．但是，我们可以使用多项式回归估计的方法．

另外，我们也可以使用支持向量机回归[9,10]的方法来解决问题1和问题2. 因为，支持向量机的方法对于具有小样本的数据估计问题也具有很好的效果．由

这样，我们就可以采用相应的方法，分别求解这两个问题． 模型假设

（1）假设对于固定的年度x 所收获的两类海产品都是按照一定的概率密度函数产生的．

（2）在未来的年度，这样的统计规律也不发生太大的变化． 模型建立

为了能够对问题的中数据变化趋势有一个清楚地直观感觉，我们将这些数据用Excel 画在坐标系中进行观察．

可以看出，我们不能用线性回归的方法来求解．下面，根据我们刚才的分析，首先采用多项式回归的方法来建模，然后再用支持向量机回归的方法来建模．

模型1 为了讨论问题的方便，我们对年度重新编号为x , x , , x ，另外，

1

2

11

给蓝鱼和蓝蟹分别编号为1，2．我们采用五次多项式回归估计（当然，可以采用其它阶数的回归多项式）．设回归函数的近似形式是如下J k 次多项式

y =∑a (j k )x j , k =1, 2 (4.4.4)

j =0J k

模型2 由于这里的数据较少，用支持向量机回归的方法是最合适的．就是要寻找一个回归函数

*k k k

=∑(αk y j -αj )K (x , x j )+b , k =1, 2

(4.4.5) k x

j =1l k

这里，l k 是第k 类海产品的样本数，K x (1), x (2)是称为核函数，其选择方法可以参考文献[9-10]，这里，我们选择径向基核函数

*k

而αk j , αj *k k αk j , αj , b

()

K x (), x (

1

(

2)

)

12

x ()-x ()

=e

-

σ (4.4.6)

(k =1,2, j =1,2, ,11)都是非负数，其意义见参考文献[10]．其中，(k =1, 2, j =1, 2, ,11)是下面优化问题的最优解

αi *∈R l , αi ∈R l

min

1l k *k

αi k *-αi k )(αk (∑j -αj )K (x i , x j )2i , j =1

l k

k *i

+ε∑(α

i =1

+α

k i

)-∑y (α

i i =1

l k

k *i

-αi k )

s .. t

∑(α

i =1

l k

(4.4.7)

k *

i

-αi k )=0

C

, k =1, 2, i , j =1, 2, , l k l k

0≤αi k *, αi k ≤

ε是事先选定的一个正数，它确定了回归函数(4.4.5)与样本函数的差别大小．详细的思想请参看文献[9,10]．我们选定的支持向量方法是解决模式识别和回归估计问题的通用方法，是建立在三大统计定律上的现代非参数统计学习方法（见文献[9,10]）．我们不需要回归函数或者识别函数的太多的信息，只要这些数据就可以了，算法会将包含在数据中的信息提取出来而用于预测或者模式识别．这种方法对于小样本问题同样适用． 模型求解

模型１的求解 我们关键是如何选择a (j k ), j =1,2, , J k , k =1,2．显然，最小二乘法的思想是一个不错的选择．建立下面的最优化问题

⎛J k (k )j (k )⎫min a t -∑ ∑j l t l ⎪, k =1,2 (4.4.8) (k )

a j , j =0,1, , J k l =1j =0

⎝⎭

根据极值的必要条件，我们得到，回归多项式满足的代数方程为

J k L k

⎛L k i +j ⎫(k )i

, t l =i 0, 1, , J k =, k 1, 2 (4.4.9) ∑a j ∑t l ⎪=∑y l

j =0l =1l =1⎝⎭

我们通过MA TLAB 编程，运行后，得到蓝鱼和蓝蟹的预测多项式分别是

y =-4.8424+6.8984x -2.4424x 2+0.4036x 3-0.0325x 4+0.0012x 5 (4.4.10)

y =-1.1729+1.7572x -0.6332x 2+0.1614x 3-0.0175x 4+0.0006x 5 (4.4.11)

L k

2

用指数函数和多项式拟合的方法，可以得到蓝鱼和蓝蟹的预测公式分别是

y =5.28571⋅(1.4635), x =1,2, ,11, （4.4.12）

y =x =1,2, ,11, . （4.4.13）

x

将原始数据与预测值分别画在同一坐标系中，可以观察到一些现象．

结果发现，用多项式预测具有随机波动的数值具有很大的偏离实际问题的本意（如对蓝蟹的多项式预测函数），在后面的时段的预测效果可能让人难以接受，即对于长期预测的效果可能比较差．但是对于短期的预测效果还是比较好的．

为此，我们可以采用用于处理预测的当前的流行方法，即基于支持向量机的回归预测的方法[9,10]．这就是我们采用模型2的原因之一．

模型2的求解 模型2涉及一个高级的模式识别和回归估计的方法[9-10]．我们直接求解优化问题，并将上述的两种方法预测的结果与支持向量机回归预测得到结果进行比较（图4.4.2和图4.4.3）．

对蓝鱼得到的预测函数（核函数中选择的σ2=5.0000002）是 y =63.60-123.28K (x , x 1)+172.68K (x , x 2)-236.90K (x , x 3)+186.13K (x , x 3)-145.70K (x , x 5)+109.91K (x , x 7)

-265.79K (x , x 8)+278.76K (x , x 9)-227.47K (x , x 10)+64.44K (x , x 11)

2

对蓝蟹得到预测函数（核函数中选择的σ=5.6）为

(4.4.14)

y =5.2083-4.9843K (x , x 1)-2.1754K (x , x 5)+0.5785K (x , x 7)-0.3656K (x , x 9)-0.8640K (x , x 10)+0.6628K (x , x 11)

(4.4.15)

结果分析

从图4.4.2和图4.4.3，我们清楚地发现，本文对于蓝蟹的多项式预测公式对原数据的拟合显然优于文献[11]的根式函数的预测结果．我们在实际问题中应该尝试使用不同次幂的回归多项式，以达到最佳的拟合．通过尝试我们发现3次多项式回归可能要更好点．同学们通过自己编程，体验研究的乐趣．

对于蓝鱼模型和蓝蟹模型，我们得到的预测结果画在图4.4.5中．

结果分析

从图4.4.4和图4.4.5可以看出，支持向量机回归的方法得到的结果最好，对于本问题来说，多项式回归预测的方法不比文献[11]的方法好．但是，如果选择合适的多项式的次数，也许会得到较好的预测，希望有兴趣的同学试试．对于没有支持向量机理论和方法的大学生来说，基于最小二乘法的多项式回归还是比较合适的．当然，有兴趣的同学可以参看文献[9-10]学习支持向量机的理论和方法． 研究性问题

读者可以尝试选择合适的回归多项式的阶使得预测更合理，或者根据观察的数据散点图，选用你认为更好的函数类型进行拟合．能否依据所给数据采用微分方程建模方法求解预测问题，或者利用最近几次历史数据值或预测值，预测以后较近时段的数据．这些都是非常有意思的问题．你会从中体会到研究的乐趣．

§4-5 随机模拟问题

对于研究对象的数量关系过于复杂或提出的解释性（定性或定量）模型难

以处理时，研究者很难得到一个能充分说明问题的符号分析模型，但又必须对研究对象的行为（随机依赖关系或者确定性关系）做出预报时，研究者可以在某种给定条件下进行多次重复的实验来收集数据，以获得这样的随机依赖关系．这种方法称为随机模拟方法[7,11]．前面对于变量之间的随机依赖关系的预报是直接利用给定数据，采用某种对回归函数的近似估计来实现的．但是，在没有这些数据情况下，我们只能采用模拟实验的方法．

在许多实际问题中，具体地进行实验来获得所需要数据是不切实际的．比如，为了确定人类对某种药物的敏感性，我们可以用小白鼠或者猴子进行模拟试验；为了能够获得人体各个器官对失重环境的适应性，我们可以进行模拟太空失重环境；为了测试电梯的某种运行方式是否合理（如停偶数上层还是停奇数层），我们不能在各种运行方式下进行多次实验，这样对顾客多有惊扰．这里的几个例子，前两个是可以有替代的试验对象，后一个则没有．在这样的情况下，我们必须设计出能够模拟实际环境或者条件的理论上的模拟仿真实验，来分析研究对象的随机依赖关系或者确定性关系．这里介绍的前两个例子也是模拟，它们是一种真实环境的模拟．而后者是借助于计算机仿真的模拟，这种模拟方法通常称为蒙特卡洛（Monete Carlo）方法．这里仅介绍这种方法．

蒙特卡洛（Monete Carlo ）方法分为确定性行为模拟和随机行为模拟．我们分别举例说明．

问题1：曲线下的面积计算－确定性问题

问题描述

我们要计算由曲线y =f (x ), (x ∈[a , b ])与直线x =a 、x =b 以及x 轴所围成

[11]

的曲边梯形的面积A ，如图4.5.1． 模型的建立

所求的面积为

⎰

b

a

f (x ) dx (4.5.1)

这个面积可以通过下面的分析给出求解的近似公式：在矩形[a , b ; Q , M ]中随机产生点P (x , y )（通过产生随机数来获得x , y ），统计出落在曲线下方的随机点

数N f 与落在整个矩形区域的随机点数N [a , b ; Q , M ]，它们的比值应当近似地等于曲边梯形的面积S f 与矩形的面积S [a , b ; P , M ]，即

S f S [a , b ; QM ]

≈

N f N [a , b ; Q , M ]

(4.5.2)

所以

S f ≈

S [a , b ; Q , M ]N f N [a , b ; Q , M ]

(4.5.3)

这里，M ≥max f (x )．

x ∈[a , b ]

问题1求解 我们采用如下的方法产生

b y =M ⋅t (4.5.4) x =a ⋅t +(1-)t ；

这里，t 为[0,1]区间均匀分布的随机数．然后，按照下面的程序进行模型的仿真试验，最后得到问题的求解． 蒙特卡洛计算面积的方法

输入 模拟中产生的随机点数N [a , b ;0, M ]． 输出 曲边梯形的面积S f 近似值． 第一步 初始化计数器counter =0；

第二步 对i =1,2, , N [a , b ; P , M ]，进行第三到第五步． 第三步 计算随机坐标(x i , y i )满足a ≤x i ≤b ;0≤y i ≤M ． 第四步 对随机坐标x i 计算f (x i )．

第五步 如果y i ≤f (x i )，则counter =counter +1．

第六步 根据公式(4.5.3)计算曲边梯形面积的近似值．

sin x

下面我们给出f (x )=在区间[0,π]上对应的曲边梯形的面积，得到的

x

蒙特卡洛近似值序列（如表4.5.1）．

sin x

表4.5.1 区间[0,π]上曲线f (x )=下的面积的蒙特卡洛近似

x

点数 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200

面积近似值 1.5708 1.7593 1.7593 1.9478 1.7090 2.0420 2.1183 1.8692 1.8361 1.8598 1.9192 1.8483 1.8173 1.9927 1.8891 1.7318 1.8923 2.0176 1.8023 1.8473 1.7443 1.6422

点数 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400

面积近似值 1.8740 1.8247 1.8297 1.7955 1.8431 1.9141 1.8026 1.8221 1.8606 1.7907 1.8240 1.7926 1.9334 1.8780 1.9019 1.8651 1.8431 1.8598 1.8543 1.8386 1.7827 1.8721

将计算结果按照迭代计算的次数递增的顺序画在图4.5.2中，明显地发现算法的

收敛性．各次模拟计算得到的近似值的平均值是1.8459，可以认为是面积的较为精确的值．通过级数的方法得到

结果分析

从结果看，蒙特卡洛方法能够较好地逼近最好的解，但是需要的计算量却是巨大的．

π

∞π∞π2k +1(-1)(-1)2k sin x

=⎰∑x dx =∑≈1.8519 (4.5.5)

0x k =02k +1! k =02k +1! 2k +1k

k

⎰

问题2：投掷一枚非均匀的骰子——随机问题

问题的描述

在概率论发展历史上非常有名的投掷硬币的试验验证了事件频率的稳定性，为概率论的发展鉴定了理论发展的基础．投掷骰子也是一个很有趣的实验． 我们的目的是计算投掷非均匀的一枚骰子，计算出现某一面的概率．总共有六个基本事件：出现第i 面就是事件A i ．用x 表示出现的点数，则事件A i 对等于随机变量x 取i 值．

我们已经知道事件的经验概率（表4.5.2）

`

表4.5.2 投掷一枚非均匀骰子出现不同面的频率 骰子出现的值

1 2 3 4 5 6

出现的频率

0.1 0.1 0.2 0.3 0.2 0.1

显然，我们可以定义事件的概率函数 f (x )=p x , x =1,2,3,4,5,6，我们的目的是利用经验频率求出这个表达式． 模型建立

显然，我们设计出独立重复的实验，用计算机来仿真完成问题的求解．我们来计算事件的频率，用事件的频率来近似代替概率．从而，给出概率函数的近似形式．即

f (i )≈p (i )=

N i

N

问题2求解 根据经验频率，适当定义下面函数以确定某个面出现的事件（表4.5.3）．

表4.5.3 出现某个面的事件与随机数x i 落在某个区间的对应关系

x i 的值

赋值

(0,0.1]

1

(0.1,0.2]

2

(0.2,0.4]

3

(0.4,0.7]

4

(0.7,0.9]

5

(0.9,1.0]

6

根据类似上面的蒙特卡洛方法，我们得到模拟概率值与经验概率值的关（表4.5.4）．

表4.5.4 蒙特卡洛方法计算事件发生的概率

骰子出现值

1 2

100

1000

5000

10000

40000

期望结果 0.1 0.1

0.1400 0.1010 0.1026 0.0992 0.0993 0.0700 0.1050 0.0958 0.0954 0.1004

3 4 5 6

0.2000 0.1940 0.1938 0.1941 0.2019 0.3700 0.3030 0.2980 0.3100 0.2952 0.1600 0.1940 0.2062 0.1989 0.2007 0.0600 0.1030 0.1036 0.1024 0.1025

0.2 0.3 0.2 0.1

结果分析

结果表明，随着模拟次数的增加，各次结果将趋向于期望的结果．

§4-6 排队服务问题

问题描述

日常生活中经常会遇到排队等候服务的现象：车站售票处排队购买车票，医院病人排队就诊，超市排队付款等等．在服务系统中，乘客、病人、顾客等统称为“顾客”，售票员、医生、收款员等统称为“服务员”．试建立模型描述服务系统的规律，并建立描述服务系统效率的数量指标，给出提高效率的措施． 问题分析

在排队系统中，服务对象的到达时间和服务时间都是随机的．研究随机服

务模型的数学工具是排队论．排队论通过对每个个别的随机服务现象的统计研究，找出反映这些随机现象平均特性的规律，从而为设计新的服务系统和改进现有服务系统的工作提供依据．

对于排队服务系统， 顾客常常注意排队的人是否太多，等候的时间是否太长，而服务员则关心他空闲的时间是否太短．于是人们常用排队的长度、等待的时间及服务利用率等指标来衡量系统的性能． 模型假设

随机服务过程包括三个部分：顾客到达规律，服务时间，排队规则．下面就这三个部分作出假设：

1. 顾客到达规律：若将时间段划分的足够细，则在[t，t+∆t ]内至多到达1个顾客，到达1个顾客的概率与时间段的长度∆t 成正比，比例系数为λ，称为到达强度，到达2个及2个以上顾客的概率为ο(∆t ) ；在彼此不相交的时间区间内到达的顾客数相互独立；

2. 服务时间：设在已经服务过一定时间的条件下再服务一段时间的概率，与一开始就服务这段时间的概率相等；

3. 排队规则：按照顾客到达的先后顺序服务，即先到先服务． 模型建立及求解

在随机服务过程中人们最关心的是队长和等待时间．由于顾客到达时刻和服务时间的随机性，队长和等待时间也是随机的．但在稳定状况下（t →∞），它们存在平稳分布．下面就队长的平稳分布，平均队长及平均等待时间建立模型．

1. 平稳分布

按照建立随机人口模型的方法考虑此问题：顾客进入服务系统相当于出生，服务完毕离开服务系统相当于死亡，这里，出生和死亡与系统人数无关．

由假设１，[0,t ]内到达的顾客数服从Poisson 分布，即有k 个顾客到达的概率为

P {X t =k }=e

-λt

(λt )

k !

k

, k =0,1, 2, ( 4.6.1)

由假使２，记Z 为服务时间，则

P {Z >τ+t Z >τ}=P {Z >t } (4.6.2) 由概率论知识知，具有这一性质的分布是指数分布．近一步假设单位时间内被服务顾客的平均人数为μ（称平均服务率），则

P {Z >t }=e -μt （4.6.3）

记服务系统时刻t 有n 个顾客（１人被服务，n -1人等待）的概率是P n (t )，则由随机人口模型得

⎧dP n (t )

=λP n -1(t )+μP n -1(t )-(λ+μ) P n (t ), n =1,2, ⎪

⎪dt

⎨ (4.6.4)

dP t ⎪0()=μP t -λP t

1()0()⎪⎩dt

在稳定状态下，P n (t ) 与t 无关，上式化为

⎧λP n -1+μP n +1-(λ+μ)P n =0⎪ ⎨

⎪⎩μP -λP 0=0

⎛λ⎫

解得P n = ⎪P 0,

⎝μ⎭

n

(4.6.5)

n =1,2, ，其中，P n 满足∑P n =1．

n -0

∞

λ

，称为服务强度，是单位时间内顾客到达数与被服务数之比．当μ

到达数小于被服务数之比时，队列长度才是有限的，故设ρ

令ρ=

所以，稳态下系统有n 个顾客的概率为

P n =(1-ρ)ρn , n =0,1,2, (4.6.6) 其中，P 0=1-ρ为服务员空闲的概率．

2. 平均队长 平均队长为

L 1=∑nP n =∑n (1-ρ)ρn =

n =0

n =0

∞

∞

ρ

1-ρ

(4.6.7)

由上式可看出，服务强度ρ越大，平均队长L 1越大，且当ρ→1, L 1→∞．

3. 平均等待时间

由假设1，顾客到达间隔服从参数为λ的指使数分布，由假设2，顾客服务时间服从参数为μ的指数分布．记顾客等待时间（含服务时间）为随机变量Y ，则Y 服从参数为λ-μ的指数分布，即

P (Y >t )=e -(μ-λ)t (4.6.8) 故平均等待时间为

W 1=与平均队长L 1相比，有

W 1=

L 1

1ρ

= (4.6.9) μ-λλ1-ρλ

(4.6.10)

即平均等待时间等于平均队长除以单位时间到达的平均人数． 研究性问题

1. 当顾客平均到达率λ增加使服务强度ρ增加，致使平均队长L 1太大时，应增加服务员．下面讨论有2个服务员的情形．

2个服务员的服务形式有两种：1个队列（M/M/2模型）；2个队列（2个M/M/1模型），下面讨论两种情形的服务效果：

对M/M/2模型

服务过程的平均服务率比原来增加1倍，为2μ，服务强度是ρ2=

λ

2μ

故由上面的方法可求出平均队长

2ρ2

L 2= (4.6.11) 2

1-ρ2平均等待时间

W 2=

L 2

=

2ρ2

λ

λ1-ρ

22

(4.6.12)

对2个M/M/1模型

λ1' λ' '

平均到达率λ1=，服务强度ρ1=与ρ2相同，故每个队列的平均长度为

μ2

ρ2'

L 1= (4.6.13)

1-ρ2

两个队列的平均长度为

' 2L 1=

2ρ2

(4.6.14) 1-ρ2

每个顾客的平均等待时间为

W =

' 1

' L 1

λ1'

=

2ρ2

(4.6.15)

λ1-ρ2与L 2, W 2相比

' 2L 1W 1'

==1+ρ2 (4.6.16) L 2W 2

由于0

W 1' 1

所以，对有多个服务员的服务系统，若从等待时间考虑，应让顾客排一个队．

对3个服务员的服务系统，可类似讨论．

2. 若从经济的角度考虑，将顾客的等待时间及服务员的服务时间视为成本，可考虑建立以此成本为目标的优化模型，以此来确定最优的服务员数目．