三、 参数估计
在确定了数据分布类型之后,下一步就是估计这个分布的参数值。目前有很多统计与数据处理软件包都能够进行常见分布的参数估计。本节列举了部分随机分布参数估计值以供参考。
为了将分布类型缩减为一个特定的随机分布并检验产生的假设,分布参数的数值估计是必需的。下表给出了在仿真中经常使用的分布的建议估计。
常用仿真随机分布参数估计量
分布 |
参数 |
估计量 |
泊松分布 |
||
指数分布 |
||
正态分布 |
||
对数正态分布 |
||
威布尔分布 |
||
伽马分布 |
||
贝塔分布 |
四、拟合优度检验
通过上面的步骤之后,已经选定了数据服从的分布类型,同时估计出该类分布的参数对应值,即确定了数据服从的概率密度函数。拟合优度检验就是评估确定后的概率密度函数表示数据输入模型的适宜程度。主要采用卡方检验和科尔莫戈洛夫-斯米尔诺夫检验两种方法,这两种方法的具体步骤参看伪随机数及随机变量产生与检验部分,下面以列举一些使用卡方检验的示例。
4.1 卡方检验泊松分布
以例2车辆到达数据为例。
步骤:
Step1:从数据的直方图来看它服从泊松分布;
Step2:通过计算样本期望,可以估计参数;
Step3:可以设计如下的假设:
H0:随机变量服从参数的泊松分布;
H1:随机变量不服从参数的泊松分布。
Step4:分组和计算见下表
其中:泊松分布的pmf如下式:
xi |
fi |
pi |
npi |
npi-fi |
(npi-fi)^2/npi |
0 |
12 |
0.01 |
1.34 |
-15.87 |
35.29 |
1 |
11 |
0.06 |
5.79 |
||
2 |
20 |
0.12 |
12.48 |
-7.52 |
4.54 |
3 |
17 |
0.18 |
17.93 |
0.93 |
0.05 |
4 |
10 |
0.19 |
19.31 |
9.31 |
4.49 |
5 |
8 |
0.17 |
16.65 |
8.65 |
4.49 |
6 |
8 |
0.12 |
11.96 |
3.96 |
1.31 |
7 |
7 |
0.07 |
7.36 |
0.36 |
0.02 |
8 |
7 |
0.04 |
3.97 |
-9.99 |
14.25 |
9 |
5 |
0.02 |
1.90 |
||
10 |
4 |
0.01 |
0.82 |
||
11 |
1 |
0.00 |
0.32 |
||
|
|
|
|
Sum |
64.44 |
其中:因为单项npi小于5,所以x=0和1进行合并;x>=8进行合并。最终分组数k=8,变量为1,则该卡方检验的自由度为6。检验统计量为64.44。
Step5:结论
使用显著性水平0.05,查卡方分布表,得临界值=12.6,远小于64.44,所以拒绝原假设,即数据不服从参数的泊松分布。因此,分析人员需要寻找一种更合适的模型或使用经验分布来进行拟合。
4.2 卡方检验均匀分布
在数的前800位小数中,数字0,1,…9出现的次数如下:
数字 |
频数 |
利用检验法,检验这些数字服从均匀分布()
解 此均匀是离散型的均匀分布,讨论的共十个数,各个数落在同一位置上的概率为0.1,共计800个位置,理论频数为=80,没有未知参数。
|
|
|
|
|
|
|
0 1 2 3 4 5 6 7 8 9 |
74 92 83 79 80 73 77 75 76 91 |
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 |
80 80 80 80 80 80 80 80 80 80 |
-6 12 3 -1 0 -7 -3 -5 -4 11 |
36 114 9 1 0 49 9 25 16 121 |
0.45 1.8 0.1125 0.0125 0 0.6125 0.1125 0.3125 0.2 1.525 |
5.125
查表
接受,服从均匀分布
4.3 卡方检验指数分布
指数分布检验步骤:
l 对样本进行升序排序,并获得最大和最小样本值[Max,Min]和样本距离D;
l 确定区间数ZoneNum;
l 获得样本分组间距GroupDistance=D/ZoneNum;
l 根据GroupDistance将排序后的样本值分为ZoneNum个组,分别为[Min,Min+GroupDistance),......,[Max-GroupDistance,Max];
l
获取每组的样本数量;
l
计算每组的理论cdf值,也就得出了
值;
l
计算,如果有些<5,将其与相邻的组进行合并再计算,获得有效组数ZoneNum2;
l
计算,计算
;
l
查表获得的值,并比较它与
的大小;如果前者大于后者,则接受随机样本来自于指数分布的假设;否则拒绝假设。
例:记录机床两次故障间隔时间如下表所示,检验机床故障间隔时间是否服从指数分布,显著水平取0.05。
8.04 |
4.77 |
12.42 |
9.22 |
3.91 |
8.50 |
14.74 |
9.20 |
5.82 |
28.56 |
10.33 |
4.62 |
7.66 |
7.24 |
39.22 |
2.14 |
15.46 |
16.16 |
16.45 |
1.18 |
22.38 |
3.12 |
0.59 |
3.76 |
3.01 |
4.00 |
4.59 |
7.07 |
1.55 |
29.58 |
10.38 |
9.59 |
12.41 |
20.05 |
0.78 |
9.34 |
8.59 |
2.54 |
1.34 |
11.50 |
1.34 |
3.82 |
9.47 |
15.83 |
14.00 |
25.77 |
3.59 |
3.62 |
0.08 |
2.87 |
21.83 |
18.41 |
21.70 |
21.95 |
1.91 |
5.97 |
2.53 |
9.21 |
2.20 |
4.86 |
0.96 |
18.79 |
6.33 |
7.16 |
5.80 |
14.65 |
0.86 |
4.33 |
4.31 |
0.07 |
0.04 |
1.33 |
17.15 |
22.81 |
6.59 |
5.94 |
3.51 |
4.22 |
5.49 |
13.27 |
26.06 |
0.29 |
1.04 |
16.29 |
1.66 |
5.81 |
3.78 |
37.63 |
0.44 |
1.77 |
2.24 |
19.26 |
9.32 |
16.39 |
11.56 |
2.67 |
8.16 |
2.40 |
3.53 |
5.14 |
解:
Step1: 计算均值为:9.018,设定假设:
H0:数据服从的指数分布;
H1:数据不服从的指数分布
Step2:设定数据分组为10,组间距离为:3.922,计算均值为:9.018,其他计算见下表。
No |
Min |
Max |
fi |
pi |
npi |
npi-fi |
(npi-fi)^2/npi |
1 |
0.000 |
3.922 |
36 |
0.353 |
35.270 |
-0.730 |
0.015 |
2 |
3.922 |
7.845 |
21 |
0.228 |
22.830 |
1.830 |
0.147 |
3 |
7.845 |
11.767 |
15 |
0.148 |
14.778 |
-0.222 |
0.003 |
4 |
11.767 |
15.690 |
7 |
0.096 |
9.566 |
2.566 |
0.688 |
5 |
15.690 |
19.612 |
9 |
0.062 |
6.192 |
-2.808 |
1.273 |
6 |
19.612 |
23.535 |
6 |
0.040 |
4.008 |
|
|
7 |
23.535 |
27.457 |
2 |
0.026 |
2.594 |
||
8 |
27.457 |
31.380 |
2 |
0.017 |
1.679 |
||
9 |
31.380 |
35.302 |
0 |
0.011 |
1.087 |
||
10 |
35.302 |
39.225 |
2 |
0.020 |
1.995 |
||
|
|
合并 |
12 |
0.114 |
11.364 |
-0.636 |
0.036 |
|
|
|
|
|
|
Sum |
2.162 |
Step3:查表获得,所以接受假设H0。
五、选择无数据的输入模型
不论希望与否,总有时候我们无法得到模型输入所需要的可靠数据,例如:实际系统根本不存在、收集数据太昂贵、数据支离破碎或者没有相应的协作关系去收集数据。面对这种情况,可能不得不依靠假设或猜想来设计输入模型。对于这些输入模型必须要在一定程度上进行输出的敏感性分析,以便把握最终结果的可信程度。
在得不到数据的情况下,一般依靠如下方法获得过程的信息:
工程数据:通常情况下,一个产品或过程都会有一个厂家提供的性能等级,例如:一个磁盘驱动器的平均故障时间是10000小时;一个激光打印机每分钟可以打印12页;一把刀具的切削速度是1厘米/秒等等。根据这些产品的过程数据,使用固定一个中心值的方法来设计输入模型。
专家建议:同在该过程或相似过程上有经验的人进行交谈。他们经常会提供一些乐观的、悲观的或最可能的时间。他们还可能告诉你这个过程是接近于常量,还是剧烈变化的,有时他们还可能定义可变性的来源。
物理或惯例的限制:大多数实际过程在性能上都有物理的限制,例如:计算机数据输入的速度受到扫描设备和传输设备速度的限制而有一个最快传输速度限制,人员在车间的行驶速度不可能超出人的正常行进速度,人工作业过程的时间受到标准操作工时的控制不能波动幅度太大等等。
对无法获得数据的过程进行输入模型设计时,可以精选一些确定性数据作为输入模型,而对于具有随机波动的过程经常需要采用随机分布作为输入模型。通常情况下,可以先考虑指数分布、三角分布、截断正态分布和均匀分布作为选择项。这几种分布的参数便于理解,并且可以反映比较广泛的数据特征,如下表所示。
可能采用的无数据输入模型
分布 |
参数 |
特征 |
可用实例 |
指数分布 |
均值 |
变化幅度大 左边有界 右边无界 |
到达间隔时间 机器无故障时间 |
三角分布 |
最小值、众数、最大值 |
对称或非对称 两边都有界 |
作业时间 |
截断正态分布 |
最小值、均值、方差、最大值 |
对称或非堆成 两边都有界 |
作业时间 |
均匀分布 |
最小值、最大值 |
所有数值等可能出现 两边都有界 |
对过程几乎不了解 |
如果时间数据的变化是独立的,并且波动比较大,则指数分布可能是不错的选择,它通常用来表示到达间隔时间,例如顾客进入餐厅、到达仓库的需求。
如果时间数据代表某项作业,而且存在一个“最可能”时间,其他的时间在其上下波动,则通常使用三角分布,因为它可以很好的反映数据小幅度或者大幅度的变化,而且它的参数很容易理解。三角分布的参数有最小值、众数和最大值,这三个参数可以很明了的估计一项活动所需时间的变化特征。三角分布的优点是允许数据在众数周围非对称分布,这种情况在实际中很普遍。三角分布也是一个有界分布:没有数据可以小于最小值,也没有数据能大于最大值。
对于可以使用三角分布表示的输入过程,也可以使用截断正态分布。截断正态分布的参数有最小值、均值、方差和最大值,数据可以对均值是对称的,也可以是非对称的,而且波动幅度可以通过方差进行控制。
为什么不使用正态分布呢?正态分布为“钟形曲线”,参数为平均值和标准差,返回的数值在平均值两边对称分布,而且是没有界的,这意味着正态分布有时会产生很小的值,有时又会产生很大的数值。当分布所代表的变量为时间数据时,正态分布返回负数将不符合实际情况,Witness会以0来代替该负值。如果正态分布的均值接近与0,则返回值有一半是0值。因此在模型中使用正态分布表示时间数据的输入模型是不合适的。如果数据集合同正态分布拟合的相当好,那么使用其他的分布以避免负值的出现,如威布尔分布、伽马分布、对数正态分布、埃尔朗分布、贝塔分布,甚至经验分布,同样能达到几乎差不多的拟合效果。
最后,如果对系统过程真的不怎么了解,但是可以猜测变量可能的最大值和最小值,可以选择均匀分布作为输入模型。
六、拟合非平稳泊松过程
非平稳泊松过程是一种比较特殊的情况,但是很值得重视,因为实际中经常会发生,是有效描述系统行为的重要方式。许多系统都依赖于外部到达,如一些服务系统、电话呼叫中心以及有外部顾客需求的制造系统等,到达过程会随着时段的不同有很大差异。例如午餐时间赶制各种汉堡包,在下午的中间时段对技术支持热线的呼叫会更加繁重,在某些季节对生产制造系统存在大量的产品需求。此时需要处理两个问题:怎样估计或确定到达率函数,以及在仿真中建立到达模式。
很多方法可以从数据中估计或确定到达率函数,这里介绍一种比较简单的方法:分段定速函数法,这种方法首先确定到达率比较平稳的时段长度,例如,一个呼叫中心的到达情况在每半小时的时段内可能比较平稳,但是在不同的时段会不一样。然后累加每个时段内的到达数量,对每个时段计算不同的到达率。
考虑三个营业员在星期一8:30~10:30的两小时内邮件到达数量,计算分段到达速率如下表。
时段 |
到达数 |
估计平均到达率 件/人/小时 |
||
People1 |
People2 |
People3 |
||
8:30-9:00 |
10 |
14 |
11 |
24 |
9:00-9:30 |
18 |
16 |
20 |
36 |
9:30-10:00 |
32 |
30 |
28 |
60 |
10:00-10:30 |
20 |
15 |
19 |
36 |
对于非平稳泊松过程在Witness中可以通过Part:Active with profile方式设置或者Distribution元素设置来实现。