# 输入建模
LLMS 索引： [llms.txt](/llms.txt)
---
输入模型为仿真提供驱动力。在排队系统的仿真中，典型的输入模型是到达间隔时间分布与服务时间分布。对于供应链仿真，输入模型包括需求与提前期的分布。对于可靠性问题的仿真，部件的失效时间分布就是输入模型的一个例子。
在真实世界的仿真应用中，为输入数据选择合适的分布，从时间和资源需求的角度看是一项主要任务。无论分析人员有多么高明，有缺陷的输入模型都会导致输出在解释时可能产生误导性的建议。
建立有用的输入数据模型包含四个步骤：
1. 从所关注的真实系统中收集数据。这往往需要投入大量的时间与资源。遗憾的是，在某些情况下无法收集数据（例如，时间极为有限、输入过程尚不存在，或者法律或规章禁止收集数据）。当数据不可得时，必须依靠专家意见和对过程的了解作出有根据的猜测。
2. 确定一个概率分布来表示该输入过程。当数据可得时，这一步通常从绘制数据的直方图（histogram）开始。给定直方图以及对过程的结构性认识，就可以选择一个分布族。所幸的是，若干为人熟知的分布在实际应用中常常能提供良好的近似。
3. 选择参数，以确定该分布族的一个具体实例。当数据可得时，这些参数可以由数据估计得出。
4. 评估所选分布及其相关参数的拟合优度（goodness of fit）。拟合优度既可以通过图形方法非正式地评估，也可以通过统计检验正式地评估。卡方检验与 Kolmogorov–Smirnov 检验是标准的拟合优度检验。若不满意所选分布对数据的近似程度，分析人员就回到第 2 步，选择另一个分布族，并重复这一过程。若该过程的若干次迭代都无法在某种假定分布形式与所收集数据之间取得拟合，则可以采用该分布的经验形式。
本章逐一讨论上述每一步。尽管现在已有大量软件可用来完成第 2、3、4 步——包括 ExpertFit® 和 \(\mathrm{S t a t} {: :} \mathrm{F i t} ^ {\circledR}\) 这样的独立程序，以及 Arena 的 Input Processor 和 @Risk 的 BestFit® 这样的集成程序——但理解软件究竟做了什么仍然重要，这样才能恰当地使用软件。遗憾的是，当两个或多个所关注的变量之间存在关系时，或者当没有数据可用时，用于输入建模的软件就不那么容易获得了。这两个主题在本章靠后的部分讨论。
## 1 数据收集
在数学、物理、化学以及其他技术类学科的教科书中，每章末尾都附有习题，供读者练习。年复一年地做这些习题，会给读者一种印象，仿佛数据是现成可得的。事实远非如此。数据收集是求解实际问题时最大的任务之一，也是仿真中最重要、最困难的问题之一。
输入建模的焦点在于把概率分布拟合到数据的统计层面，这些分布随后将为仿真提供驱动输入。然而，在"拟合"发生之前，需要有准确且相关的输入数据可用（或可收集），并且需要充分理解数据的性质。人们常常把数据直接交给输入建模软件，而软件总会给数据拟合出某个概率模型，不论这种拟合是否有正当理由。下面的例子说明了这种做法常见的陷阱：
### 例 1：陈旧数据
一项仿真研究将考察如何缩短患者在癌症筛查诊所中耗费的时间。筛查流程的一部分由护士采集患者病史。在 2002 年的一项护理活动研究中，一家咨询公司收集了护士采集患者病史的观测数据。这家医院仍保存着这些数据。当前项目是否应当用它们来做输入建模？
在这种情况下，答案很可能是否定的，因为当前项目的工程师们深入调查后发现，2008 年护士已停止在纸质病历上记录筛查笔记，改为直接把内容录入笔记本电脑。因此，2002 年的数据对应的是另一种不同的病史采集过程，可能会严重歪曲当前的过程。如果由于所需时间或费用而无法收集当前数据，那么宁可使用护士本人给出的估计，也不要依赖那份过时的研究。第 6 节描述了把专家意见转化为概率分布的方法。工程师们还可以主动让自己的病史被采集一次并计时，用得到的几个实际数据点来检验护士的估计。
### 例 2：出乎意料的数据
出于安保目的，进入某政府大楼的人员必须通过一道包含金属探测器的安检。在实施新的安保流程之前，将开展一项仿真研究以评估其影响。为便于研究，收集了 1000 个通过金属探测器站点所需时间（以分钟计）的观测值。这些数据的均值约为 1/2 分钟，标准差几乎与之相同。这提示通过金属探测器的时间可以用指数分布作为良好的输入模型，因为指数分布的均值与标准差总是相等，而且在排队系统中指数分布常被用来建模服务时间。
然而，分析人员对这一输入模型感到不安，因为根据她的经验，通过金属探测器只需几秒钟。她绘制了数据的直方图，如图 1 所示，这立刻表明通过金属探测的时间由两个明显不同的分布构成。仔细思考后，她意识到有些人毫无意外地通过探测器，而另一些人则触发了警报，必须取出身上残留的金属物品再试一次，某些情况下甚至需要接受人身检查。
**图 1 以分钟计量的 1000 个通过金属探测器时间观测值的直方图**
因此，这批数据需要三个不同的输入模型：一个是未触发探测器而通过所需时间的模型（直方图中最左侧的 764 个观测值）；一个是警报确实响起时通过所需时间的模型（直方图中最右侧的 236 个观测值）；还有一个是一个人触发警报的概率模型。虽然指数分布当然能与数据的汇总统计量（均值与标准差）相匹配，但拟合出的指数分布会遗漏安保检查的一个关键特征，几乎必然导致不准确的仿真。
遗憾的是，仔细检查数据后，分析人员又发现了一个复杂情况：图 1 的直方图显示样本中存在负值（共九个）。通过金属探测的时间为负显然是不可能的，那么分析人员应当怎么做？可能是观测记录本身正确，只是在计算时间时出了错；如果是这样，分析人员或许能够重建正确的数值。若不能，则在拟合分布之前应当舍弃这九个负值。但在估计某人触发警报的概率时，是否也应当舍弃它们呢？也就是说，估计值应当是 236/1000 还是 236/991？如果分析人员能够说服自己，这九个负值实际上是记录错误的小时间值，那么她可能希望在估计该概率时保留它们。如果她不确定，那么仿真应当分别用这两个估计各运行一次，看看是否存在实质性的影响。
### 例 3：时变数据
一家计算机制造商设有客户支持中心，产品所有者可以打电话咨询问题或提出疑虑。该呼叫中心的坐席服务时间从东部时间上午 8:00 到晚上 8:00。由于坐席的工资与福利是主要成本，因此计划开展一项仿真研究以帮助确定人员配置水平。呼叫中心拥有详尽的每日呼叫量记录，显示他们平均每天接到 12,277 通电话。固定时间或空间内随机事件的数量常常可以用泊松分布很好地建模，于是对数据拟合了一个泊松分布（使用本章后面描述的方法），并且看起来拟合良好。当到达数服从泊松分布时，到达之间的时间间隔——这正是仿真呼叫所需要的——必定服从指数分布。该指数分布到达率的自然估计是每 12 小时工作日 12,277 次到达，但由于"一天"是相当大的时间单位，它被换算为每分钟 17 通电话。
上一段的分析看起来是可靠的，泊松分布的选择（意味着到达间隔时间服从指数分布）甚至得到了数据的支持。然而，数据的聚合层级——一个 12 小时工作日——掩盖了呼叫中心一个非常重要的特征：呼叫的到达率可能而且通常确实会随时间（一天中的时段、一周中的星期几，甚至一年中的时节，例如假日送礼季刚结束时）而显著变化。我们常常假定所关注的过程是平稳的（stationary），即其基本性质不随时间变化。然而，如果存在趋势，那么考虑这些趋势可能远比选择用哪个概率分布来表示系统中剩余的不确定性更为重要。在第 5 节中，我们描述一种拟合非平稳泊松过程的方法，它很可能更适合这个呼叫中心。
### 例 4：相关数据
某大型零售商已从订单数据库中提取出某个品牌运动饮料的月度需求，用于一项供应链仿真。针对这批数据的一个简单输入模型是：月度需求独立同分布（i.i.d.）。本章描述了对此类数据拟合分布并评估拟合效果的方法。但同样重要的是，还要检验"独立同分布"这一假设本身。
**图 2 带时间滞后的需求数据**
检验这些数据是否同分布的一种非正式方法是按时间顺序绘制需求，以查找趋势（例如需求随时间下降）或季节性（存在规律的需求周期，例如每年夏季需求高、随后冬季需求低）。带有趋势或季节性的数据通常应当先"去趋势"（通常通过对趋势或季节模式拟合一个模型），然后再对剩余变异拟合一个分布。
独立性可以通过带时间滞后的散点图来检验。图 2 是一幅散点图，把该运动饮料每个月的 demand 与下一个月的 demand 对应作图（例如 1 月的需求对 2 月的需求）。显然存在很强的负相关，这意味着某个月需求大，下个月需求往往就小，反之亦然。这一现象可能是零售商某个月过量备货、随后在下个月少订货以消化多余库存所导致的。大于一个月的时间滞后同样可以考察。如果盲目地采用独立同分布输入这一标准假设，就会漏掉这一重要现象，从而导致供应链管理效果不佳。在第 7 节中，我们给出对相关输入数据建模的方法。
正如这些例子所示，仅仅拥有数据并不足以支撑有效的输入建模。数据可能过时，也可能"脏"（含有错误）。有时，把数据转换成可用形式或"清洗"掉错误所需的工作量或成本，可能与获取数据本身相当。输入建模几乎总是要求分析人员既运用自己的判断，也使用恰当的统计工具。了解哪些输入模型最可靠、哪些最不可靠，对于判断仿真研究所得结论的可靠性非常重要。而且，由于关于输入模型正确性的不确定性永远无法彻底消除，明智的做法是用若干个看似合理的输入模型分别运行仿真，看看结论是稳健的，还是对这些选择高度敏感。
从一次真实的数据收集经历中可以学到许多教训。本章末尾的前五道习题给出了一些情境，学生可以在其中获得这类经验。
下面的建议可能有助于加强和便利数据收集工作，尽管它们并不包罗万象。
1. 把时间花在规划上是值得的。这可以从一次练习或预观测开始。试着在预观测时收集数据。为此设计表格。这些表格很可能在正式收集数据之前必须修改多次。留意异常情形，并考虑如何处理它们。在可能的情况下，对系统录像，事后通过观看录像提取数据。即使数据将自动收集（例如通过计算机收集数据），规划也很重要，以确保获得恰当的数据。当数据已由他人收集时，务必留出充足的时间把数据转换成可用的格式。
2. 尽量在收集数据的同时分析数据。弄清楚正在收集的数据是否足以提供仿真所需的分布。查明正在收集的某些数据是否对仿真无用。没有必要收集多余的数据。
3. 尽量合并同质的数据集。检查相继时间段内以及相继日期的同一时间段内数据的同质性。例如，检查下午 2:00 至 3:00 与下午 3:00 至 4:00 数据的同质性，并检查星期四与星期五下午 2:00 至 3:00 的数据是否同质。在检查同质性时，初步检验可以看分布的均值（例如平均到达间隔时间）是否相同。两样本 t 检验可用于此目的。更彻底的分析要求检验分布是否等价，也许可以通过分位数-分位数图（quantile-quantile plot，后文描述）来进行。
4. 注意数据删失（censoring）的可能性，即所关注的量未被完整观测。这一问题最常出现在分析人员关注完成某个过程所需时间时（例如生产一个零件、治疗一位患者，或某个部件失效），但该过程在观测期开始之前就已开始，或在观测期结束之后才结束。删失会导致特别长的加工时间被排除在数据样本之外。
5. 要发现两个变量之间是否存在关系，可绘制散点图。有时只要看一看散点图，就能判明所关注的变量之间是否存在关系。
6. 考虑这样一种可能：一组看起来独立的观测序列实际上具有自相关（autocorrelation）。自相关可以存在于相继的时间段中，也可以存在于相继的顾客之间。
7. 牢记输入数据与输出或性能数据之间的区别，并确保收集的是输入数据。输入数据通常表示那些很大程度上超出系统控制范围的不确定性量，不会因为为改进系统所作的变更而改变。另一方面，输出数据表示系统在
受到这些输入作用时的性能表现，而这可能正是我们试图改进的东西。在排队仿真中，顾客到达时间通常是输入，而顾客延误则是输出。性能数据对模型确认很有用。
再次说明，这些只是几点建议。一般而言，必须以极为审慎的态度对待数据收集与分析。
## 2 用数据识别分布
在本节中，我们讨论当数据可得并且被认为是独立同分布时，选择输入分布族的方法。分布族中具体分布的确定通过估计其参数来完成，如第 3 节所述。第 6 节讨论数据不可得的情形。
### 2.1 直方图
频数分布或直方图（histogram）有助于识别分布的形状。直方图按如下方式构造：
1. 把数据的取值范围划分为若干区间。（区间通常等宽；不过，如果对频数的高度作相应调整，也可以使用不等宽的区间。）
2. 标注横轴，使之与所选区间一致。
3. 找出每个区间内出现的频数。
4. 标注纵轴，以便为每个区间绘制总出现次数。
5. 在纵轴上绘制频数。
类别区间（class interval）的数目取决于观测值的个数以及数据中散布或离散的程度。Hines 等人 [2002] 指出，在实践中，把类别区间的数目取为大致等于样本量的平方根通常效果不错。如果区间过宽，直方图就会很粗糙或呈块状，其形状和其他细节就表现不好。如果区间过窄，直方图就会参差不齐，无法把数据平滑化。同一批数据的参差、粗糙与恰当的直方图示例见图 3。现代数据分析软件通常允许方便、交互地改变区间大小，直到找到合适的选择为止。
连续数据的直方图对应于理论分布的概率密度函数。如果数据是连续的，则穿过每个类别区间频数中心点所画出的线应当形成类似 pdf 的形状。
对于数据点很多的离散数据，其直方图应当为数据取值范围内的每一个值设一个格子。然而，如果数据点很少，可能就有必要合并相邻格子，以消除直方图的参差外观。如果直方图对应的是离散数据，它看起来应当像一个概率质量函数。
**图 3 参差、粗糙与恰当的直方图：(a) 原始数据——过于参差；(b) 合并相邻格子——过于粗糙；(c) 合并相邻格子——恰当**
### 例 5：离散数据
在 20 周期间的五个工作日中，对上午 7:00 至 7:05 这 5 分钟内到达某交叉口西北角的车辆数进行了监测。表 1 给出了所得数据。表中第一项表示有 12 个 5 分钟时段内零辆车到达，有 10 个时段内一辆车到达，依此类推。
**表 1 5 分钟时段内的到达数**
| 每时段到达数 | 频数 | 每时段到达数 | 频数 |
|---|---|---|---|
| 0 | 12 | 6 | 7 |
| 1 | 10 | 7 | 5 |
| 2 | 19 | 8 | 5 |
| 3 | 17 | 9 | 3 |
| 4 | 10 | 10 | 3 |
| 5 | 8 | 11 | 1 |
**图 4 每时段到达数的直方图**
汽车数量是离散变量，而且数据充足，因此直方图可以为数据取值范围内的每一个可能值设一个格子。所得直方图如图 4 所示。
### 例 6：连续数据
对随机抽取的电子元件样本在 1.5 倍额定电压下进行了寿命试验，并记录其寿命（或失效时间），以天为单位：
| 79.919 | 3.081 | 0.062 | 1.961 | 5.845 |
|---|---|---|---|---|
| 3.027 | 6.505 | 0.021 | 0.013 | 0.123 |
| 6.769 | 59.899 | 1.192 | 34.760 | 5.009 |
| 18.387 | 0.141 | 43.565 | 24.420 | 0.433 |
| 144.695 | 2.663 | 17.967 | 0.091 | 9.003 |
| 0.941 | 0.878 | 3.371 | 2.157 | 7.579 |
| 0.624 | 5.380 | 3.148 | 7.078 | 23.960 |
| 0.590 | 1.928 | 0.300 | 0.002 | 0.543 |
| 7.004 | 31.764 | 1.005 | 1.147 | 0.219 |
| 3.217 | 14.382 | 1.008 | 2.336 | 4.562 |
寿命通常被视为连续变量，这里记录到小数点后三位。把数据放入类别区间即可绘制直方图。数据的取值范围相当大，从 0.002 天到 144.695 天。然而，大多数数值（50 个中的 30 个）落在 0 至 5 天的范围内。取区间宽度为 3，得到表 2。然后用表 2 的数据绘制出图 5 所示的直方图。
**图 5 元件寿命的直方图**
### 2.2 选择分布族
绘制直方图的目的是推断出某个已知的 pdf 或 pmf。一个分布族的选择，依据是所考察情境中可能出现的分布以及直方图的形状。因此，如果已经收集了到达间隔时间数据，且直方图的形状类似于指数分布，就可以假定为指数分布。类似地，如果正在测量货盘上货物的重量，而直方图看起来关于均值对称，那么假定为正态分布就是有充分理由的。
**表 2 电子元件数据**
| 元件寿命（天） | 频数 |
|---|---|
| \(0 \leq x_{j} < 3\) | 23 |
| \(3 \leq x_{j} < 6\) | 10 |
| \(6 \leq x_{j} < 9\) | 5 |
| \(9 \leq x_{j} < 12\) | 1 |
| \(12 \leq x_{j} < 15\) | 1 |
| \(15 \leq x_{j} < 18\) | 2 |
| \(18 \leq x_{j} < 21\) | 0 |
| \(21 \leq x_{j} < 24\) | 1 |
| \(24 \leq x_{j} < 27\) | 1 |
| \(27 \leq x_{j} < 30\) | 0 |
| \(30 \leq x_{j} < 33\) | 1 |
| \(33 \leq x_{j} < 36\) | 1 |
| . | . |
| . | . |
| . | . |
| \(42 \leq x_{j} < 45\) | 1 |
| . | . |
| . | . |
| . | . |
| \(57 \leq x_{j} < 60\) | 1 |
| . | . |
| . | . |
| . | . |
| \(78 \leq x_{j} < 81\) | 1 |
| . | . |
| . | . |
| . | . |
| \(144 \leq x_{j} < 147\) | 1 |
指数分布、正态分布和泊松分布是经常遇到的，从计算角度看也不难分析。虽然 beta 分布、伽马分布和 Weibull 分布分析起来较为困难，但它们提供的形状极为丰富，在对底层概率过程建模时不应被忽视。也许起初假定的是指数分布，却发现它对数据拟合不佳。下一步应当是考察拟合不佳出现在何处。如果拟合不佳出现在分布的某一尾部，那么伽马分布或 Weibull 分布或许能更充分地拟合数据。
人们已经构造出的概率分布实际上有成百上千种；其中许多是带着某种特定的物理过程而创建的。选择分布的一个辅助手段是利用分布的物理基础作为指引。下面是一些例子：
二项分布（Binomial）：建模 n 次试验中成功的次数，其中各次试验相互独立且具有共同的成功概率 p；例如，在一批 n 个芯片中发现的有缺陷芯片的数量。
负二项分布（Negative Binomial，包含几何分布）：建模为取得 k 次成功所需的试验次数；例如，为找到 4 个有缺陷芯片而必须检查的芯片数量。
泊松分布（Poisson）：建模在固定时间或空间内发生的独立事件的个数；例如，1 小时内到达某商店的顾客数，或在 30 平方米金属薄板上发现的缺陷数。
正态分布（Normal）：建模一个可视为若干分量过程之和的过程的分布；例如，装配一件产品的时间就是各项装配作业所需时间的总和。注意正态分布允许取负值，而对于加工时间来说负值可能是不可能的。
对数正态分布（Lognormal）：建模一个可视为若干分量过程之积（即相乘）的过程的分布；例如，当利息按复利计算时，一项投资的收益率就是若干期收益的乘积。
指数分布（Exponential）：建模独立事件之间的时间间隔，或具有无记忆性（memoryless）的加工时间（知道已经过去多少时间，并不能提供关于还要多久过程才完成的任何信息）；例如，来自大量相互独立行动的潜在顾客的到达之间的间隔时间。指数分布是变异性很高的分布；它有时被过度使用，因为它常常导致数学上易于处理的模型。回顾一下，如果事件之间的时间间隔服从指数分布，那么固定时间段内事件发生的个数就服从泊松分布。
伽马分布（Gamma）：一种极为灵活的分布，用于建模非负随机变量。伽马分布可以通过加上一个常数而从 0 处平移开。
Beta 分布（Beta）：一种极为灵活的分布，用于建模有界（具有固定的上限和下限）的随机变量。Beta 分布可以通过加上一个常数而从 0 处平移开，也可以通过乘以一个常数而获得比 [0, 1] 更大的取值范围。
厄兰分布（Erlang）：建模可以看作若干个指数分布过程之和的过程；例如，当一台计算机和两台备份计算机都失效时计算机网络失效，而每台计算机的失效时间都服从指数分布。厄兰分布是伽马分布的一个特例。
Weibull 分布（Weibull）：建模部件的失效时间；例如磁盘驱动器的失效时间。指数分布是 Weibull 分布的一个特例。
离散或连续均匀分布（Uniform）：建模完全的不确定性：所有结果等可能。在没有数据时，这一分布常常被不恰当地使用。
三角分布（Triangular）：建模只知道分布的最小值、最可能值和最大值的过程；例如测试一件产品所需的最小、最可能和最大时间。这一模型通常比均匀分布有显著改进。
经验分布（Empirical）：从实际收集到的数据中重抽样；通常在似乎没有合适的理论分布时使用。
在选择分布时，不要忽视过程的物理特征。这个过程天然是离散取值还是连续取值？它是有界的还是没有天然界限？这类不依赖数据的知识，有助于缩小可供选择的分布族的范围。还要记住，对任何随机输入过程而言，并不存在"真实"的分布。输入模型是对现实的一种近似，因此目标是从仿真实验中获得一个能产生有用结果的近似。
建议读者完成习题 6 至习题 11，以进一步了解本节提到的各种分布的形状。观察参数变化时形状的改变，是很有启发性的。
### 2.3 分位数-分位数图
如第 2.1 节所讨论的直方图的构造，以及如第 2.2 节所讨论的分布形状的识别，都是为一批样本数据选择分布族的必要组成部分。然而，直方图在评估所选分布的 \(\mathscr {f} t\) 方面并不那么有用。当数据点个数较少（比如 30 个或更少）时，直方图可能相当参差。此外，我们对拟合程度的感知取决于直方图区间的宽度。而且，即便区间选得很好，把数据分组到格子中也使得很难把直方图与连续概率密度函数作比较。分位数-分位数图（quantile-quantile plot，\(( q \ - q )\) 图）是评估分布拟合的有用工具，它不受这些问题的困扰。
如果 X 是 cdf 为 F 的随机变量，那么 X 的 q 分位数（quantile）就是满足 \(F ( \gamma ) =\)\(P ( X \leq \gamma ) = q\) 的值 \(\gamma\)，其中 \(0 < q < 1\)。当 F 存在逆函数时，我们写作 \(\gamma = F^{- 1} ( q )\)。现在令 \(\{x_{i} , i = 1 , 2 , \dots , n \}\) 是来自 X 的一个数据样本。把观测值从小到大排序，并记为 \(\{y_{j} , j = 1 , 2 , . . . n \}\)，其中 \(y_{1} \leq y_{2} \leq \cdot \cdot \cdot \leq y_{n}\)。令 j 表示秩或序号。因此，最小者 \(j = 1\)，最大者 \(j = n\)。\(q {-} q\) 图基于这样一个事实：\(y_{j}\) 是 X 的 \(( j - 1 / 2 ) / n\) 分位数的一个估计。换言之，
$$
y_{j} \text{is approximately} F^{- 1} \left(\frac {j - \frac {1}{2}}{n}\right)
$$
现在，假设我们选择了 cdf 为 \(F\) 的分布作为 X 分布的一种可能表示。如果 \(F\) 是某个恰当分布族的成员，那么把 \(y_{j}\) 对 \(F^{- 1} ( ( j - 1 / 2 ) / n )\)) 作图将近似为一条直线。如果 F 来自恰当的分布族且参数值也恰当，那么该直线的斜率将为 1。另一方面，如果所假定的分布不恰当，这些点就会偏离直线，
通常是以一种系统性的方式偏离。是否拒绝某个假设模型的判断是主观的。
### 例 7：正态 Q-Q 图
一台机器人用于在装配线上安装汽车车门。人们认为安装时间服从正态分布。该机器人能够精确测量安装时间。机器人自动采集了 20 个安装时间的样本，结果如下，数值以秒为单位：
| 99.79 | 99.56 | 100.17 | 100.33 |
|---|---|---|---|
| 100.26 | 100.41 | 99.98 | 99.83 |
| 100.23 | 100.27 | 100.02 | 100.47 |
| 99.55 | 99.62 | 99.65 | 99.82 |
| 99.96 | 99.90 | 100.06 | 99.85 |
样本均值为 99.99 秒，样本方差为 (0.2832)^2 秒^2。这些值可以作为正态分布的均值与方差的参数估计。现将观测值从小到大排序如下：
| j | 值 | j | 值 | j | 值 | j | 值 |
|---|---|---|---|---|---|---|---|
| 1 | 99.55 | 6 | 99.82 | 11 | 99.98 | 16 | 100.26 |
| 2 | 99.56 | 7 | 99.83 | 12 | 100.02 | 17 | 100.27 |
| 3 | 99.62 | 8 | 99.85 | 13 | 100.06 | 18 | 100.33 |
| 4 | 99.65 | 9 | 99.90 | 14 | 100.17 | 19 | 100.41 |
| 5 | 99.79 | 10 | 99.96 | 15 | 100.23 | 20 | 100.47 |
然后把排序后的观测值对 \(F^{- 1} ( ( j - 1 / 2 ) / 20 )\) 作图，其中 \(j = 1 , 2 , \dots , 20\)，F 是均值为 99.99、方差为 (0.2832)^2 的正态分布的 cdf，从而得到 q-q 图。所绘制的点连同数据的直方图一起显示在图 6 中。注意，从直方图很难判断数据是否能被正态分布很好地表示，但在 q-q 图中呈一条直线的总体观感却相当清晰，并支持正态分布的假设。
在评估 q-q 图的线性程度时，应当考虑以下几点：
a. 观测值永远不会精确地落在一条直线上。
b. 排序后的值并不独立；它们已被排过秩。因此，如果一个点位于某条直线上方，那么下一个点很可能也位于该直线上方。而且，这些点不太可能均匀地散布在直线两侧。
c. 极端值（最大值和最小值）的方差远高于图中部值的方差。在极端处可以接受较大的偏差。图中部诸点的线性比极端处的线性更重要。
**图 6 安装时间的 q-q 图与直方图**
现代数据分析软件通常包含生成 q-q 图的工具，尤其是针对正态分布的。\(q {-} q\) 图也可用于比较两个数据样本，看它们能否用同一个分布表示（即它们是否同质）。\(\operatorname{I f} x_{1} , x_{2} , \ldots , x_{n}\) 是随机变量 \(X ,\) 的一个样本，而 \(z_{1} , z_{2} , \ldots , z_{n}\) 是随机变量 \(Z ,\) 的一个样本，那么把 X 的排序值与 Z 的排序值对应作图，如果两个样本都能被同一个分布很好地表示，则图形将近似呈一条直线（Chambers、Cleveland 与 Tukey [1983]）。
## 3 参数估计
选定了分布族之后，下一步就是估计该分布的参数。本节描述许多有用分布的估计量。此外，现在已有许多软件包——其中一些已集成到仿真语言中——可以用来计算这些估计。
### 3.1 预备统计量：样本均值与样本方差
在许多情况下，样本均值，或样本均值与样本方差，被用来估计所假设分布的参数；参见例 7。下面的段落给出三组用于计算样本均值和样本方差的公式。当可得的是离散或连续的原始数据（raw data）时，使用式 (1) 和式 (2)。当数据是离散的且已被归入频数分布时，使用式 (3) 和式 (4)。当数据是离散或连续的且已被放入类别区间时，使用式 (5) 和式 (6)。式 (5) 和式 (6) 是近似式，只有在原始数据不可得时才应使用。
如果容量为 n 的样本中的观测值为 \(X_{1} , X_{2} , \ldots , X_{n}\)，则样本均值 \(\bar {X}\) 定义为
$$
\bar {X} = \frac {\sum_{i = 1} ^ {n} X_{i}}{n}\tag{1}
$$
样本方差 \(S^{2}\) 定义为
$$
S^{2} = \frac {\sum_{i = 1} ^ {n} X_{i} ^ {2} - n \bar {X} ^ {2}}{n - 1}\tag{2}
$$
如果数据是离散的且已被归入频数分布，则可以对式 (1) 和式 (2) 加以修改，从而大大提高计算效率。样本均值可以计算为
$$
\bar {X} = \frac {\sum_{j = 1} ^ {k} f_{j} X_{j}}{n}\tag{3}
$$
样本方差为
$$
S^{2} = \frac {\sum_{j = 1} ^ {k} f_{j} X_{j} ^ {2} - n \bar {X} ^ {2}}{n - 1}\tag{4}
$$
其中 k 是 X 的不同取值的个数，\(f_{j}\) 是 X 的取值 \(X_{j}\) 的观测频数。
### 例 8：分组数据
对表 1 中的数据进行分析可得 \(n = 100 , f_{1} = 12 , X_{1} = 0 , f_{2} = 10 , X_{2} = 1 , \ldots ,\)\(\textstyle \sum_{j = 1} ^ {k} f_{j} X_{j} = 364\)，以及 \(\begin{array} {r} {\sum_{j = 1} ^ {k} = f_{j} X_{j} ^ {2} = 2080} \end{array}\)。由式 (3) 得，
$$
\bar {X} = \frac {364}{100} = 3. 64
$$
并由式 (4) 得，
$$
S^{2} = \frac {2080 - 100 (3 . 64) ^ {2}}{99} = 7. 63
$$
样本标准差 S 就是样本方差的平方根。在本例中，\({\boldsymbol {S}} =\)\(\sqrt {7 . 63} = 2 . 76\)。用式 (1) 和式 (2) 也会得到关于 \(\bar {X}\) 和 \(S^{2}\) 完全相同的结果。当取值是连续的时，如果可能，最好使用原始数据。然而，有时收到的数据已经被放进了类别区间。这时就不再可能得到精确的样本均值和样本方差。在这种情况下，样本均值和样本方差由下列公式近似：
$$
\bar {X} \doteq \frac {\sum_{j = 1} ^ {c} f_{j} m_{j}}{n}\tag{5}
$$
以及
$$
S^{2} \doteq \frac {\sum_{j = 1} ^ {c} f_{j} m_{j} ^ {2} - n \bar {X} ^ {2}}{n - 1}\tag{6}
$$
其中 \(f_{j}\) 是第 j 个类别区间中的观测频数，\(m_{j}\) 是第 j 个区间的中点，c 是类别区间的个数。
### 例 9：类别区间中的连续数据
假设例 6 中元件寿命的原始数据或者已被丢弃，或者已经丢失。然而，表 2 中的数据仍然可得。为了近似 X¯ 和 \(S^{2}\) 的值，使用式 (5) 和式 (6)。创建如下数值：\(f_{1} = 23 , m_{1} = 1 . 5 , f_{2} = 10 \ / \mathrm{\Omega}\)\(m_{2} = 4 . 5 , \ldots , \sum_{j = 1} ^ {49} f_{j} m_{j} = 614\) 以及 \(\textstyle \sum_{j = 1} ^ {49} f_{j} m_{j} ^ {2} = 37 , 226 . 5\)。取 \(n = 50\)，由式 (5) 近似得 X 为
$$
\bar {X} \doteq \frac {614}{50} = 12. 28
$$
然后，由式 (6) 近似得 \(S^{2}\) 为
$$
S^{2} \doteq \frac {37 , 226 . 5 - 50 (12 . 28) ^ {2}}{49} = 605. 849
$$
以及
$$
S \doteq 24. 614
$$
对例 6 中的原始数据应用式 (1) 和式 (2)，得到 \(\bar {X} \ = \ 11 . 894\) 和 \(S = 24 . 953\)。因此，当原始数据被丢弃或丢失时，可能会产生误差。
**表 3 仿真中常用分布的建议估计量**
| 分布 | 参数 | 建议估计量 |
|---|---|---|
| 泊松 | \(\alpha\) | \(\widehat{\alpha} = \bar{X}\) |
| 指数 | \(\lambda\) | \(\widehat{\lambda} = \frac{1}{\bar{X}}\) |
| 伽马 | \(\beta, \theta\) | \(\widehat{\beta}\)（见表 A.9） \(\widehat{\theta} = \frac{1}{\bar{X}}\) |
| 正态 | \(\mu, \sigma^{2}\) | \(\widehat{\mu} = \bar{X}\)
\(\widehat{\sigma}^{2} = S^{2}\)（无偏） |
| 对数正态 | \(\mu, \sigma^{2}\) | \(\widehat{\mu} = \bar{X}\)（对数据进行 ln 变换后） \(\widehat{\sigma}^{2} = S^{2}\)（对数据进行 ln 变换后） |
| Weibull（\(\nu = 0\)） | \(\alpha, \beta\) | \(\widehat{\beta}_{0} = \frac{\bar{X}}{S}\)
\(\widehat{\beta}_{j} = \widehat{\beta}_{j-1} - \frac{f(\widehat{\beta}_{j-1})}{f'(\widehat{\beta}_{j-1})}\) 见式 (11) 与式 (14) 中 \(f(\widehat{\beta})\) 与 \(f'(\widehat{\beta})\) 的定义 迭代至收敛 \(\widehat{\alpha} = \left( \frac{1}{n} \sum_{i=1}^{n} X_{i}^{\widehat{\beta}} \right)^{1/\widehat{\beta}}\) |
| Beta | \(\beta_{1}, \beta_{2}\) | \(\Psi(\widehat{\beta}_{1}) + \Psi(\widehat{\beta}_{1} - \widehat{\beta}_{2}) = \ln(G_{1})\)
\(\Psi(\widehat{\beta}_{2}) + \Psi(\widehat{\beta}_{1} - \widehat{\beta}_{2}) = \ln(G_{2})\) 其中 \(\Psi\) 为双伽马函数（digamma function），\(G_{1} = \left( \prod_{i=1}^{n} X_{i} \right)^{1/n}\) 且 \(G_{2} = \left( \prod_{i=1}^{n} (1 - X_{i}) \right)^{1/n}\) |
### 3.2 建议的估计量
需要分布参数的数值估计，以便把分布族缩小为一个具体的分布，并检验由此产生的假设。表 3 给出了仿真中常用分布的建议估计量。除了对正态分布的 \(\sigma^{2}\) 估计作了去除偏倚的修正之外，这些估计量都是基于原始数据的最大似然估计量。（如果数据位于类别区间中，这些估计量必须加以修改。）关于均匀分布、二项分布和负二项分布的参数估计，请读者参阅 Fishman [1973] 和 Law [2007]。三角分布通常在无数据可用时使用，其参数由对最小值、最可能值和最大可能值的有根据猜测获得；如果只有最小值和最大值可用，均匀分布也可以这样使用。
下面各段给出这些估计量的使用示例。读者应当牢记，参数是一个未知的常数，而估计量是一个统计量（或随机变量），因为它依赖于样本值。为了在这里清楚地区分二者，如果某个参数记作 α，则其估计量将记作 α。
### 例 10：指数分布
我们将用指数分布来说明获得参数估计量的最大似然法。速率参数为 λ 的指数分布的密度函数为
$$
f (x) = \lambda e^{- \lambda x}
$$
其中 \(x \ge 0\)。如果 \(X_{1} , X_{2} , \ldots , X_{n}\) 独立同分布，且各自服从该指数分布，那么它们的联合分布为
$$
f (x_{1}, x_{2}, \ldots , x_{n}) = \lambda e ^ {- \lambda x_{1}} \lambda e ^ {- \lambda x_{2}} \dots \lambda e ^ {- \lambda x_{n}} = \lambda^{n} e ^ {- \lambda \sum_{i = 1} ^ {n} x_{i}}\tag{7}
$$
联合分布大致表示任意一组 \(x_{i}\) 取值的似然。
最大似然法把这个关系反转过来：如果我们有一个数据样本，比如 \(X_{1} , X_{2} , \ldots , X_{n}\)，但 λ 未知，那么最大似然估计量（MLE）就是使该样本似然（由式 (7) 度量）最大化的那个 λ 值。因此，指数分布 λ 的 MLE 使如下似然函数最大化
$$
L (\lambda) = \lambda^{n} e ^ {- \lambda \sum_{i = 1} ^ {n} X_{i}}
$$
求 MLE 通常是一项关于未知参数最大化似然函数的微积分练习。由于似然函数是密度函数的乘积且取值为正，有时最大化似然函数的对数会更容易，在本例中即为
$$
\ln L (\lambda) = n \ln (\lambda) - \lambda \sum_{i = 1} ^ {n} X_{i}
$$
对 \(\lambda\) 求导并令其等于 \(_ 0\)，得到
$$
\frac {n}{\lambda} - \sum_{i = 1} ^ {n} X_{i} = 0
$$
于是得到 MLE \(\begin{array} {r} {\widehat {\lambda} = n / \sum_{i = 1} ^ {n} X_{i} = 1 / \bar {X}} \end{array}\)，其中 \(\bar {X}\) 是观测值的样本均值。
例如，假设例 6 中的数据来自指数分布，则由 \(\bar {X}\) 得到 λ 的 MLE 为
$$
\widehat {\lambda} = \frac {1}{\bar {X}} = \frac {1}{11 . 894} = 0. 084 \text{per day}
$$
### 例 11：Weibull 分布
假设已经抽取了一个容量为 n 的随机样本 \(. X_{1} , X_{2} , \ldots , X_{n} . ,\)，并且假定这些观测值来自 Weibull 分布。可以证明其似然函数为
$$
L (\alpha , \beta) = \frac {\beta^{n}}{\alpha^{\beta n}} \left[ \prod_{i = 1} ^ {n} X_{i} ^ {(\beta - 1)} \right] \exp \left[ - \sum_{i = 1} ^ {n} \left(\frac {X_{i}}{\alpha}\right) ^ {\beta} \right]\tag{8}
$$
最大似然估计是使 \(L ( \alpha , \beta )\) 最大、等价地使 ln \(L ( \alpha , \beta )\)（记作 \(l ( \alpha , \beta )\)）最大的那些 \(\widehat {\alpha}\) 和 \(\widehat {\beta}\) 值。\(l ( \alpha , \beta )\) 的最大值通过对 \(\partial l ( \alpha , \beta ) / \partial \alpha\) 和\(\partial l ( \alpha , \beta ) / \partial \beta\) 取偏导、令其各自为零、并求解所得方程而获得，代入后这些方程变为
$$
f (\beta) = 0\tag{9}
$$
以及
$$
\alpha = \left(\frac {1}{n} \sum_{i = 1} ^ {n} X_{i} ^ {\beta}\right) ^ {1 / \beta}\tag{10}
$$
其中
$$
f (\beta) = \frac {n}{\beta} + \sum_{i = 1} ^ {n} \ln X_{i} - \frac {n \sum_{i = 1} ^ {n} X_{i} ^ {\beta} \ln X_{i}}{\sum_{i = 1} ^ {n} X_{i} ^ {\beta}}\tag{11}
$$
最大似然估计\(\widehat {\alpha}\) 和 \({\widehat {\beta}} , ,\) 是式 (9) 与式 (10) 的解。首先，通过下面说明的迭代过程求得\(\widehat {\beta}\)。然后在式 (10) 中取 \(\beta = {\widehat {\beta}}\) 求得 \(\widehat {\alpha}\)。式 (9) 是非线性的，因此必须使用数值分析技术来求解。在表 \(^ {3 ,}\) 中，给出的计算 \(\widehat {\beta}\) 的建议迭代方法为
$$
\widehat {\beta} _ {j} = \widehat {\beta} _ {j - 1} - \frac {f (\widehat {\beta} _ {j - 1})}{f^{\prime} (\widehat {\beta} _ {j - 1})}\tag{12}
$$
式 (12) 在逼近 \({\widehat {\beta}} ,\) 时采用牛顿法（Newton's method），其中 \({\widehat {\beta}} _ {j}\) 是第 j 次迭代，从表 3 给出的初始估计 \({\widehat {\beta}} _ {0}\) 开始，如下：
$$
\widehat {\beta} _ {0} = \frac {\bar {X}}{S}\tag{13}
$$
如果初始估计\({\widehat {\beta}} _ {0} . ,\) 与解 \({\widehat {\beta}} , ,\) 足够接近，那么当 \(j \longrightarrow \infty\) 时 \({\widehat {\beta}} _ {j}\) 趋于 \(\widehat \beta\)。在牛顿法中，\(\widehat {\beta}\) 是通过大小为 \(f ( \widehat {\beta} _ {j - 1} ) / f^{\prime} ( \widehat {\beta} _ {j - 1} )\) 的增量逐步逼近的。用式 (11) 计算 \(f ( \widehat {\beta} _ {j - 1} )\)，用式 (14) 计算 \(f^{\prime} ( \widehat {\beta} _ {j - 1} )\)，如下：
$$
f^{\prime} (\beta) = - \frac {n}{\beta^{2}} - \frac {n \sum_{i = 1} ^ {n} X_{i} ^ {\beta} (\ln X_{i}) ^ {2}}{\sum_{i = 1} ^ {n} X_{i} ^ {\beta}} + \frac {n \left(\sum_{i = 1} ^ {n} X_{i} ^ {\beta} \ln X_{i}\right) ^ {2}}{\left(\sum_{i = 1} ^ {n} X_{i} ^ {\beta}\right) ^ {2}}\tag{14}
$$
式 (14) 可由式 (11) 对 \(\beta .\) 求导 \(f ( \beta )\) 得到。迭代过程一直持续到 \(f ( \widehat {\beta} _ {j} ) \doteq 0\)，即直到 $| f ( \widehat {\beta} _ {j} ) | \leq 0 . 001 $ 。
考虑例 6 中给出的数据。这些数据涉及电子元件的失效，看起来来自指数分布。在例 10 中，参数 \(\widehat {\lambda}\) 是在数据来自指数分布的假设下估计的。如果这一假设被拒绝，另一个假设就是数据来自 Weibull 分布。之所以怀疑是 Weibull 分布，是因为这些数据涉及电子元件失效，而失效是突然发生的。
用式 (13) 计算 \(\widehat {\beta} _ {0}\)。对于例 \(6 , n = 50 , \bar {X} = 11 . 894 , \bar {X} ^ {2} =\)\(141 . 467\) 中的数据，以及 \(\begin{array} {r} {\sum_{i = 1} ^ {50} X_{i} ^ {2} = 37 , 575 . 850 ;} \end{array}\)；于是由式 (2) 求得 \(S^{2}\) 为
$$
S^{2} = \frac {37 , 578 . 850 - 50 (141 . 467)}{49} = 622. 650
$$
且 \(S = 24 . 953\)。因此，
$$
\widehat {\beta} _ {0} = \frac {11 . 894}{24 . 953} = 0. 477
$$
要用式 (12) 计算 \({\widehat {\beta}} _ {1}\)，需要根据式 (11) 和式 (14) 计算 \(f ( \widehat {\beta_{0}} ) \mathrm{a n d} f^{\prime} ( \widehat {\beta_{0}} )\)。还需要下列附加数值：\(\textstyle \sum_{i = 1} ^ {50} X_{i} ^ {\widehat {\beta_{0}}} = 115 . 125 , \textstyle \sum_{i = 1} ^ {50} \ln X_{i} = 38 . 294\)\(\begin{array} {r} {\sum_{i = 1} ^ {50} X_{i} ^ {\widehat {\beta_{0}}} \ln X_{i} = 292 . 629} \end{array}\)，以及 \(\textstyle \sum_{i = 1} ^ {50} X_{i} ^ {\widehat {\beta_{0}}} ( \ln X_{i} ) ^ {2} = 1057 . 781\)。于是，
$$
f (\widehat {\beta_{0}}) = \frac {50}{0 . 477} + 38. 294 - \frac {50 (292 . 629)}{115 . 125} = 16. 024
$$
以及
$$
f^{\prime} \left(\widehat {\beta_{0}}\right) = \frac {- 50}{(0 . 477) ^ {2}} - \frac {50 (1057 . 781)}{115 . 125} + \frac {50 (292 . 629) ^ {2}}{(115 . 125) ^ {2}} = - 356. 110
$$
然后，由式 (12) 得，
$$
\widehat {\beta} _ {1} = 0. 477 - \frac {16 . 024}{- 356 . 110} = 0. 522
$$
四次迭代之后，$| f ( \widehat {\beta} _ {3} ) | \le 0 . 001 $，此时 \(\widehat {\beta} \doteq \widehat {\beta} _ {4} = 0 . 525 ,\) 就是式 (9) 的近似 解。表 4 列出了完成每次迭代所需的数值。
现在，可以在式 (10) 中取 \(\beta = \widehat {\beta} = 0 . 525\) 计算\(\widehat {\alpha}\)，如下：
$$
\widehat {\alpha} = \left[ \frac {130 . 608}{50} \right] ^ {1 / 0. 525} = 6. 227
$$
**表 4 Weibull 分布参数的迭代估计**
| j | \(\widehat{\beta}_{j}\) | \(\sum_{i=1}^{50} X_{i}^{\widehat{\beta}_{j}}\) | \(\sum_{i=1}^{50} X_{i}^{\widehat{\beta}_{j}} \ln X_{i}\) | \(\sum_{i=1}^{50} X_{i}^{\widehat{\beta}_{j}} (\ln X_{i})^{2}\) | \(f(\widehat{\beta}_{j})\) | \(f'(\widehat{\beta}_{j})\) | \(\widehat{\beta}_{j+1}\) |
|---|---|---|---|---|---|---|---|
| 0 | 0.477 | 115.125 | 292.629 | 1057.781 | 16.024 | -356.110 | 0.522 |
| 1 | 0.522 | 129.489 | 344.713 | 1254.111 | 1.008 | -313.540 | 0.525 |
| 2 | 0.525 | 130.603 | 348.769 | 1269.547 | 0.004 | -310.853 | 0.525 |
| 3 | 0.525 | 130.608 | 348.786 | 1269.614 | 0.000 | -310.841 | 0.525 |
如果 \({\widehat {\beta}} _ {0}\) 与 \({\widehat {\beta}} ,\) 足够接近，该过程收敛很快，通常四到五次迭代即可。 然而，如果该过程看似发散，则应对 \({\widehat {\beta}} _ {0}\) 尝试其他初值；例如，取初始估计的一半或两倍。
手工估计 Weibull 分布参数的艰巨任务，凸显了为输入建模提供软件支持的价值。
### 例 12：泊松分布
假设需要分析表 1 中的到达数据。考察图 4 提示可以假定为参数 α 未知的泊松分布。由表 3 可知，α 的估计量是 \({\bar {X}} ,\)，这已在例 8 中求出。于是，\(\widehat {\alpha} = 3 . 64\)。回顾一下，对于泊松分布，真实的均值与方差相等。在例 8 中，样本方差估计为 \(S^{2} = 7 . 63\)。然而，绝不应期望样本均值与样本方差精确相等，因为二者都是随机变量。
### 例 13：对数正态分布
某投资组合中 10 项投资的收益率分别为 18.8、27.9、21.0、6.1、37.4、5.0、22.9、1.0、3.1 和 8.3 个百分点。为了估计这些数据对数正态模型的参数，我们首先对数据取自然对数，得到 2.9、3.3、3.0、1.8、3.6、1.6、3.1、0、1.1 和 2.1。然后令 \(\widehat {\mu} = \bar {X} = 2 . 3\) 以及 \(\widehat {\sigma} ^ {2} = S^{2} = 1 . 3\)。
### 例 14：正态分布
正态分布的参数\(\mu\) 与 \(\sigma^{2}\)，由 X¯ 与 \(S^{2}\) 估计，如表 3 所示。例 7 中的 \(q {-} q\) 图得出的分布假设是安装时间服从正态分布。由式 (1) 和式 (2)，例 7 的数据给出 \(\widehat {\mu} = \bar {X} = 99 . 9865\) 和 \(\widehat {\sigma} ^ {2} = S^{2} = ( 0 . 2832 ) ^ {2}\) 秒^2。
### 例 15：伽马分布
伽马分布的估计量 \(\widehat {\beta}\) 通过查阅表 A.9 来选定。表 A.9 需要 计算量 \(1 / M\)，其中
$$
M = \ln {\bar {X}} - \frac {1}{n} \sum_{i = 1} ^ {n} \ln X_{i}\tag{15}
$$
另外，由表 3 可见 \(\widehat {\theta}\) 由下式给出
$$
\widehat {\theta} = \frac {1}{\bar {X}}\tag{16}
$$
提前期常常服从伽马分布。假设与 20 个订单相关的提前期（以天计）已被精确测量如下：
| 订单 | 提前期（天） | 订单 | 提前期（天） |
|---|---|---|---|
| 1 | 70.292 | 11 | 30.215 |
| 2 | 10.107 | 12 | 17.137 |
| 3 | 48.386 | 13 | 44.024 |
| 4 | 20.480 | 14 | 10.552 |
| 5 | 13.053 | 15 | 37.298 |
| 6 | 25.292 | 16 | 16.314 |
| 7 | 14.713 | 17 | 28.073 |
| 8 | 39.166 | 18 | 39.019 |
| 9 | 17.421 | 19 | 32.330 |
| 10 | 13.905 | 20 | 36.547 |
要估计 \(\widehat {\beta}\) 和 \({\widehat {\theta}} ,\)，首先需要由式 (15) 计算 M。这里由式 (1) 求得 \(\bar {X}\) 为
$$
\bar {X} = \frac {564 . 32}{20} = 28. 22
$$
然后
$$
\ln \bar {X} = 3. 34
$$
接着
$$
\sum_{i = 1} ^ {20} \ln X_{i} = 63. 99
$$
于是
$$
M = 3. 34 - \frac {63 . 99}{20} = 0. 14
$$
以及
$$
1 / M = 7. 14
$$
通过在附录表 A.9 中插值，得\({\widehat {\beta}} = 3 . 728 . .\)。最后，由式 (16) 得
$$
\widehat {\theta} = \frac {1}{28 . 22} = 0. 035
$$
### 例 16：Beta 分布
每个月携带商店优惠券的顾客百分比必定介于 0 与 100% 之间。在某商店八个月的观测给出如下数值：25%、74%、20%、32%、81%、47%、31% 和 8%。为了对这些数据拟合 Beta 分布，我们首先需要把所有数值除以 100，将其重新标定到区间 [0, 1] 上，得到 0.25、0.74、0.20、0.32、0.81、0.47、0.31、0.08。
参数 \(\beta_{1} , \beta_{2}\) 的最大似然估计量是表 3 中所示方程组的解。这类方程可以用现代符号/数值计算程序（如 Maple）求解；图 7 给出了一个求解 Beta 参数的 Maple 程序。在本例中，解为 \(\widehat {\beta} _ {1} = 1 . 47\) 和 \(\widehat {\beta} _ {2} = 2 . 16\)。
```txt
betaMLE := proc(X, n)
local G1, G2, beta1, beta2, eqns, solns;
G1 := product(X[i], i=1..n)^(1/n);
G2 := product(1-X[i], i=1..n)^(1/n);
eqns := {Psi(beta1) - Psi(beta1 + beta2) = ln(G1),
Psi(beta2) - Psi(beta1 + beta2) = ln(G2)};
solns := fsolve(eqns, {beta1=0..infinity, beta2=0..infinity});
RETURN(solns);
end;
```
**图 7 计算 Beta 分布参数最大似然估计的 Maple 程序**
## 4 拟合优度检验
在本节中，将 Kolmogorov–Smirnov 检验与卡方检验应用于关于输入数据分布形式的假设。
拟合优度检验为评估某个候选输入模型的适用性提供了有益的指导；然而，在实际应用中并不存在唯一正确的分布，因此你不应当成为这类检验结论的奴隶。理解样本量的影响尤其重要。如果可用数据很少，那么拟合优度检验不太可能拒绝任何候选分布；而如果可用数据很多，那么拟合优度检验很可能拒绝所有候选分布。因此，未能拒绝某个候选分布应当被视为支持该选择的一条证据，而拒绝某个输入模型也只是反对该选择的一条证据。
### 4.1 卡方检验
检验"容量为 n 的随机变量 X 的随机样本服从某种特定分布形式"这一假设的一种程序是卡方拟合优度检验。该检验把"把数据的直方图与候选密度函数或质量函数的形状作比较"这一直观想法形式化。当样本量较大、且参数由最大似然法估计时，该检验对离散与连续两种分布假设都有效。检验程序首先把 n 个观测值排列成 k 个类别区间或格子。检验统计量由下式给出
$$
\chi_{0} ^ {2} = \sum_{i = 1} ^ {k} \frac {(O_{i} - E_{i}) ^ {2}}{E_{i}}\tag{17}
$$
其中 \(O_{i}\) 是第 i 个类别区间中的观测频数，\(E_{i}\) 是该类别区间中的期望频数。每个类别区间的期望频数计算为 \(E_{i} = n p_{i}\)，其中 \(p_{i}\) 是与第 i 个类别区间相关联的理论的、假设的概率。
可以证明，\(\chi_{0} ^ {2}\) 近似服从自由度为 \(k - s - 1\) 的卡方分布，其中 s 表示由样本统计量估计的假设分布参数的个数。假设如下：\(H_{0} {\mathrm{:}}\) 随机变量 X 服从由参数估计值所给出的分布假设。\(H_{1} {\mathrm{:}}\)：随机变量 X 不服从该分布。
临界值 \(\chi_{\alpha , k - s - 1} ^ {2}\) 可在表 A.6 中查到。若 \(\chi_{0} ^ {2} > \chi_{\alpha , k - s - 1} ^ {2} .\)，则拒绝原假设 \(H_{0}\)。在应用该检验时，如果期望频数过小，\(\chi_{0} ^ {2}\) 不仅会反映观测频数对期望频数的偏离，也会反映期望频数本身过小这一点。关于 \(E_{i} ,\) 的最小取值虽然没有一致意见，但 \(3 , 4 .\), 和 5 这些值被广泛采用。如果某个 \(E_{i}\) 值过小，可以把它与相邻类别区间中的期望频数合并。相应的 \(O_{i}\) 值也应当合并，并且每合并一个格子，k 就减少 1。
如果被检验的分布是离散的，那么随机变量的每一个取值都应当是一个类别区间，除非为了满足最小期望格子频数的要求而必须合并相邻的类别区间。对于离散情形，如果不需要合并相邻格子，则
$$
p_{i} = p (x_{i}) = P (X = x_{i})
$$
否则，\(p_{i}\) 通过对适当的相邻格子的概率求和得到。
**表 5 连续数据类别区间数目的建议**
| 样本量 n | 类别区间数 k |
|---|---|
| 20 | 不要使用卡方检验 |
| 50 | 5 到 10 |
| 100 | 10 到 20 |
| >100 | \(\sqrt{n}\) 到 \(n/5\) |
如果被检验的分布是连续的，则类别区间由 \([ a_{i - 1} , a_{i} )\) 给出，其中 \(a_{i - 1}\) 和 \(a_{i}\) 是第 i 个类别区间的端点。对于假定了 pdff(x) 或假定了 cdf \(F ( x ) , p_{i}\) 的连续情形，可以计算为
$$
p_{i} = \int_ {a_{i - 1}} ^ {a_{i}} f (x) d x = F (a_{i}) - F (a_{i - 1})
$$
对于离散情形，类别区间的数目由按需合并相邻格子之后得到的格子数决定。然而，对于连续情形，类别区间的数目必须指定。虽然没有可遵循的一般规则，但表 5 中的建议有助于确定连续数据的类别区间数目。
### 例 17：对泊松假设应用卡方检验
在例 12 中，分析了例 5 给出的车辆到达数据。数据的直方图如图 \(^ {4 ,}\) 所示，看起来服从泊松分布；因此求得参数 \(\widehat {\alpha} = 3 .\).64。于是形成如下假设：\(H_{0} {\mathrm{:}}\) 随机变量服从泊松分布。\(H_{1}\)：随机变量不服从泊松分布。
泊松分布的 pmf 为：
$$
p (x) = \left\{\begin{array}{l l} \frac {e^{- \alpha} \alpha^{x}}{x !}, & x = 0, 1, \ldots \\ 0, & \text{otherwise} \end{array} \right.\tag{18}
$$
对于 \(\widehat {\alpha} = 3 . 64\)，由式 (18) 得到与 x 的各个取值相对应的概率：
$$
\begin{array}{l l} \mathrm{p(0)=0.026} & \mathrm{p(6)=0.085} \\ \mathrm{p(1)=0.096} & \mathrm{p(7)=0.044} \\ \mathrm{p(2)=0.174} & \mathrm{p(8)=0.020} \\ \mathrm{p(3)=0.211} & \mathrm{p(9)=0.008} \\ \mathrm{p(4)=0.192} & \mathrm{p(10)=0.003} \\ \mathrm{p(5)=0.140} & \mathrm{p(\geq 11) = 0.001} \end{array}
$$
**表 6 例 17 的卡方拟合优度检验**
| \(x_i\) | 观测频数 \(O_i\) | 期望频数 \(E_i\) | \((O_i - E_i)^2\)\(E_i\) | |
|---|---|---|---|---|
| 0 | 12 | 22 | 2.6 | 7.87 |
| 1 | 10 | 9.6 | | |
| 2 | 19 | 17.4 | 0.15 | |
| 3 | 17 | 21.1 | 0.80 | |
| 4 | 10 | 19.2 | 4.41 | |
| 5 | 8 | 14.0 | 2.57 | |
| 6 | 7 | 8.5 | 0.26 | |
| 7 | 5 | 17 | 4.4 | 11.62 |
| 8 | 5 | 2.0 | | |
| 9 | 3 | 0.8 | | |
| 10 | 3 | 0.3 | | |
| ≥11 | 1 | 0.1 | | |
| | 100 | 100.0 | 27.68 | |
根据这些信息构造出表 6。\(E_{1}\) 的值由 \(n p_{0} = 100 ( 0 . 026 ) =\)2.6 给出。以类似方式计算其余的 \(E_{i}\) 值。由于 \(E_{1} = 2 . 6 < 5 , E_{1}\) 与 \(E_{2}\) 被合并。在这种情况下，\(O_{1}\) 与 \(O_{2}\) 也被合并，k 减少 1。出于同样的原因，最后五个类别区间也被合并，k 再减少 4。
计算得到的 \(\chi_{0} ^ {2}\) 为 27.68。\(\chi^{2}\) 表值的自由度为 \(k - s - 1 =\)\(7 - 1 - 1 = 5\)。这里，\(s = 1\)，因为有一个参数\({\widehat {\alpha}} , ,\) 是由数据估计的。在 0.05 显著性水平下，临界值 \(\chi_{0 . 05 , 5} ^ {2}\)为 11.1。因此，\(H_{0}\) 将在 0.05 显著性水平下被拒绝。于是，分析人员可能希望寻找拟合更好的模型，或者使用数据的经验分布。
### 4.2 等概率的卡方检验
如果被检验的是连续分布假设，则应当使用概率相等而非区间宽度相等的类别区间。这已为许多作者所推荐 [Mann 与 Wald, 1942；Gumbel, 1943；Law, 2007；Stuart、Ord 与 Arnold, 1998]。应当指出，该程序不适用于以类别区间形式收集、且原始数据已被丢弃或丢失的数据。
遗憾的是，并没有一种简便的方法来确定与每个区间相关联的概率，以使给定规模的检验的功效（power）最大化。（检验的功效定义为拒绝一个错误假设的概率。）然而，如果使用等概率，则 \(p_{i} = 1 / k\)。我们建议
$$
E_{i} = n p_{i} \geq 5
$$
于是，代入 \(p_{i}\) 得到
$$
\frac {n}{k} \geq 5
$$
并对 k 求解得
$$
k \leq \frac {n}{5}\tag{19}
$$
式 (19) 被用于得出表 5 中关于最大类别区间数的建议。
如果所假设的分布是正态、指数或 Weibull 分布，本节描述的方法是直截了当的。例 18 说明了如何对指数分布实施该程序。如果所假设的分布是伽马分布（但不是厄兰分布）或某些其他分布，那么类别区间端点的计算就很复杂，可能需要对密度函数作数值积分。在这种情况下，统计分析软件非常有用。
### 例 18：指数分布的卡方检验
在例 10 中，分析了例 6 给出的失效数据。数据的直方图如图 5 所示，看起来服从指数分布，于是计算出参数 \(\widehat {\lambda} = 1 / \bar {X} = 0 . 084\)。因此形成如下假设：\(H_{0} {\mathrm{:}}\) 随机变量服从指数分布。\(H_{1} ;\)：随机变量不服从指数分布。
为了用等概率区间进行卡方检验，必须求出类别区间的端点。式 (19) 表明区间数应当小于或等于 \(n / 5 .\)。这里，\(n = 50\)，因此 \(k \leq 10\)。在表 5 中，建议使用 5 到 10 个类别区间。令 \(k = 8 ;\)，则每个区间具有概率 \(p = 0 . 125\)。每个区间的端点由指数分布的 cdf 计算如下：
$$
F (a_{i}) = 1 - e ^ {- \lambda a_{i}}\tag{20}
$$
其中 \(a_{i}\) 表示第 i 个区间的端点，\(i = 1 , 2 , \ldots , k\)。由于 \(F ( a_{i} )\) 是从零到 \(a_{i} , F ( a_{i} ) = i p\) 的累积面积，因此式 (20) 可写为
$$
i p = 1 - e ^ {- \lambda a_{i}}
$$
或
$$
e ^ {- \lambda a_{i}} = 1 - i p
$$
两边取对数并对 \(a_{i}\) 求解，得到指数分布 k 个等概率区间端点的一般结果：
$$
a_{i} = - \frac {1}{\lambda} \ln (1 - i p), \quad i = 0, 1, \dots , k\tag{21}
$$
**表 7 例 18 的卡方拟合优度检验**
| 类别区间 | 观测频数 \(O_i\) | 期望频数 \(E_i\) | \((O_i - E_i)^2\)\(E_i\) |
|---|---|---|---|
| [0, 1.590) | 19 | 6.25 | 26.01 |
| [1.590, 3.425) | 10 | 6.25 | 2.25 |
| [3.425, 5.595) | 3 | 6.25 | 0.81 |
| [5.595, 8.252) | 6 | 6.25 | 0.01 |
| [8.252, 11.677) | 1 | 6.25 | 4.41 |
| [11.677, 16.503) | 1 | 6.25 | 4.41 |
| [16.503, 24.755) | 4 | 6.25 | 0.81 |
| [24.755,∞) | 6 | 6.25 | 0.01 |
| | 50 | 50 | 39.6 |
无论 \(\lambda ,\) 取何值，式 (21) 总会有 \(a_{0} = 0\) 和 \(a_{k} = \infty\)。取 \(\widehat {\lambda} = 0 . 084\) 和 \(k = 8 , a_{1}\) 由式 (21) 计算为
$$
a_{1} = - \frac {1}{0 . 084} \ln (1 - 0. 125) = 1. 590
$$
对 \(i = 2 , 3 , \dots , 7\) 继续应用式 (21)，得到 \(a_{2} , \ldots , a_{7}\) 分别为 3.425、5.595、8.252、11.677、16.503 和 24.755。由于 \(k = 8 , a_{8} = \infty\)。第一个区间是 [0, 1.590)，第二个区间是 [1.590, 3.425)，依此类推。预期每个区间中将落入 0.125 比例的观测值。观测值、期望值以及对 \(\chi_{0} ^ {2}\) 计算值的贡献如表 7 所示。\(\chi_{0} ^ {2}\) 的计算值为 39.6。自由度由 \(k - s - 1 = 8 - 1 - 1 = 6 . \mathrm{A t} \alpha = 0 . 05\) 给出，\(\chi_{0 . 05 , 6} ^ {2}\) 的表值为 12.6。由于 \(\chi_{0} ^ {2} > \chi_{0 . 05 , 6} ^ {2} ,\)，原假设被拒绝。（\(\chi_{0 . 01 , 6} ^ {2}\) 的值为 16.8，因此在显著性水平 \(\alpha = 0 . 01\) 下原假设也会被拒绝。）
### 4.3 Kolmogorov–Smirnov 拟合优度检验
卡方拟合优度检验可以容纳由数据估计参数的情况，其结果是自由度相应减少（每估计一个参数减少一个）。卡方检验要求把数据放入类别区间；在连续分布假设的情形下，这种分组是人为的。改变类别数和区间宽度会影响计算所得和查表所得的卡方值。一种假设可能在数据按某种方式分组时被接受，而在按另一种方式分组时被拒绝。此外，卡方检验统计量的分布只是近似已知，而且该检验的功效有时相当低。基于这些考虑，人们希望有卡方检验之外的拟合优度检验。Kolmogorov–Smirnov 检验把考察 q-q 图背后的思想形式化了。
Kolmogorov–Smirnov 检验可以用来检验数值的均匀性。这两种用途都属于拟合优度检验的范畴。任何连续分布假设都可以进行拟合优度检验。
Kolmogorov–Smirnov 检验在样本量较小、且没有参数由数据估计时特别有用。当已经作出参数估计时，表 A.8 中的临界值是有偏的；具体地说，它们过于保守。在这里，"保守"意味着临界值会过大，从而导致比所设定的更小的第 I 类（α）错误。在某些情况下可以算出 α 的精确值，本节末尾将对此加以讨论。
当假定为指数分布时，Kolmogorov–Smirnov 检验不需要任何特殊的表。下面的例子说明在这种情况下如何应用该检验。（注意，在本例中不需要估计分布的参数，因此我们可以使用表 A.8。）
### 例 19：指数分布的 Kolmogorov–Smirnov 检验
假设在 100 分钟的时间区间内收集了 50 个到达间隔时间（以分钟计）（按发生顺序排列）：
| 0.44 | 0.53 | 2.04 | 2.74 | 2.00 | 0.30 | 2.54 | 0.52 | 2.02 | 1.89 | 1.53 | 0.21 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2.80 | 0.04 | 1.35 | 8.32 | 2.34 | 1.95 | 0.10 | 1.42 | 0.46 | 0.07 | 1.09 | 0.76 |
| 5.55 | 3.93 | 1.07 | 2.26 | 2.88 | 0.67 | 1.12 | 0.26 | 4.57 | 5.37 | 0.12 | 3.19 |
| 1.63 | 1.46 | 1.08 | 2.06 | 0.85 | 0.83 | 2.44 | 1.02 | 2.24 | 2.11 | 3.15 | 2.90 |
| 6.58 | 0.64 | | | | | | | | | | |
原假设及其对立假设构造如下：\(H_{0} {\mathrm{:}}\) 到达间隔时间服从指数分布。\(H_{1}\)：到达间隔时间不服从指数分布。
数据是在从 0 到 \(T = 100\) 分钟这一区间内收集的。可以证明，如果到达间隔时间 \(\{T_{1} , T_{2} , \dots \}\) 的底层分布是指数分布，那么到达时间在区间 [0, T] 上服从均匀分布。到达时间 \(T_{1} , T_{1} + T_{2} , T_{1} + T_{2} + T_{3} , . . . ,\)\(T_{1} + \cdots + T_{50}\) 是通过累加到达间隔时间得到的。然后把到达时间归一化到 [0, 1] 区间，以便应用 Kolmogorov–Smirnov 检验。在 [0, 1] 区间上，这些点将是 \([ T_{1} / T , ( T_{1} + T_{2} ) / T , \dots , ( T_{1} + \cdot \cdot \cdot + T_{50} ) / T ]\)。所得的 50 个数据点如下：
| 0.0044 | 0.0097 | 0.0301 | 0.0575 | 0.0775 | 0.0805 | 0.1059 | 0.1111 | 0.1313 | 0.1502 |
|---|---|---|---|---|---|---|---|---|---|
| 0.1655 | 0.1676 | 0.1956 | 0.1960 | 0.2095 | 0.2927 | 0.3161 | 0.3356 | 0.3366 | 0.3508 |
| 0.3553 | 0.3561 | 0.3670 | 0.3746 | 0.4300 | 0.4694 | 0.4796 | 0.5027 | 0.5315 | 0.5382 |
| 0.5494 | 0.5520 | 0.5977 | 0.6514 | 0.6526 | 0.6845 | 0.7008 | 0.7154 | 0.7262 | 0.7468 |
| 0.7553 | 0.7636 | 0.7880 | 0.7982 | 0.8206 | 0.8417 | 0.8732 | 0.9022 | 0.9680 | 0.9744 |
按照该程序得到 \(D^{+}\) 为 0.1054，\(D^{-}\) 为 0.0080。因此，Kolmogorov–Smirnov 统计量为 \(D \ = \ \operatorname* {m a x} ( 0 . 1054 , 0 . 0080 ) \ = \ 0 . 1054\)。由表 A.8 得到的
在显著性水平 \(\alpha \ : = \ : 0 . 05\) 且 \(n = 50\) 时 D 的临界值为 \(D_{0 . 05} \ =\)\(1 . 36 / \sqrt {n} = 0 . 1923\)；然而，\(D = 0 . 1054\)，因此"到达间隔时间服从指数分布"这一假设不能被拒绝。
Kolmogorov–Smirnov 检验已被改造，以便在若干种参数由数据估计的情形下使用。检验统计量的计算是相同的，但使用不同的临界值表。不同的分布假设需要不同的临界值表。Lilliefors [1967] 提出了一种正态性检验。原假设声称总体属于正态分布族中的某一个，而不指定分布的参数。有兴趣的读者不妨研读 Lilliefors 的原始工作；他描述了如何用仿真来确定临界值。
Lilliefors [1969] 还针对指数分布修正了 Kolmogorov–Smirnov 检验的临界值。Lilliefors 再次使用随机抽样来获得近似的临界值，但 Durbin [1975] 随后得到了精确分布。Connover [1998] 给出了正态分布与指数分布的 Kolmogorov–Smirnov 检验示例。他还提到了若干其他读者可能感兴趣的 Kolmogorov–Smirnov 型检验。
Anderson-Darling 检验是一种在精神上与 Kolmogorov–Smirnov 检验相似的检验。与 Kolmogorov–Smirnov 检验一样，Anderson-Darling 检验基于经验 cdf 与拟合 cdf 之间的差异；与 Kolmogorov–Smirnov 检验不同的是，Anderson-Darling 检验基于一种更综合的差异度量（不只是最大差异），并且对分布尾部的差异更为敏感。Anderson-Darling 检验的临界值同样取决于候选分布以及是否估计了参数。幸运的是，该检验和 Kolmogorov–Smirnov 检验都已在许多支持仿真输入建模的软件包中实现。
### 4.4 p 值与"最佳拟合"
要应用拟合优度检验，必须选择显著性水平。回顾一下，显著性水平是错误地拒绝 \(H_{0} {\mathrm{:}}\)：随机变量服从该分布假设的概率。传统的显著性水平是 0.1、0.05 和 0.01。在高速计算出现之前，拥有一小组标准值使得制出有用的临界值表成为可能。如今，大多数统计软件按需计算临界值，而不是把它们存进表中。因此，分析人员可以采用不同的显著性水平，比如 0.07。
然而，许多软件包并不要求预先指定显著性水平，而是为检验统计量计算一个 \(p \mathrm{.}\)-值。\(p \mathrm{-}\)-值是指对于给定的检验统计量取值，恰好会拒绝 \(H_{0}\) 的显著性水平。因此，大的 p 值往往表明拟合良好（为了拒绝，我们不得不接受很大的出错机会），而小的 p 值则表明拟合不佳（要接受，我们就必须要求几乎没有风险）。
回顾例 17，其中用卡方检验来检验车辆到达数据的泊松假设。检验统计量的值为 \(\chi_{0} ^ {2} = 27 . 68\)，自由度为 5。该检验统计量的 \(p \mathrm{-}\) 值为 0.00004，这意味着我们会在 0.00004 的显著性水平上拒绝"数据服从泊松分布"的假设。（回顾一下，我们曾在 0.05 水平上拒绝该假设；现在我们知道，即使在更低的水平上也不得不拒绝它。）
p 值可以看作拟合程度的一种度量，值越大越好。这提示我们可以把手中所有的分布都拟合一遍，为每次拟合计算一个检验统计量，然后选择产生最大 p 值的那个分布。据我们所知，并没有输入建模软件实现这一特定算法，但许多此类软件确实包含"最佳拟合"选项，即软件在评估所有可行模型之后向用户推荐一个输入模型。软件也可能考虑其他因素——例如数据是离散的还是连续的、有界的还是无界的——但最终仍会用某种拟合的汇总度量（如 p 值）来对分布排序。这样做并无不妥，但有几点需要牢记：
1. 软件可能对数据的物理基础一无所知，而这类信息却能提示哪些分布族是恰当的。（参见第 2.2 节中的清单。）记住，输入建模的目标常常是填补空白或平滑数据，而不是找到一个与给定样本尽可能贴合的输入模型。
2. 回顾一下，厄兰分布和指数分布都是伽马分布的特例，而指数分布也是更灵活的 Weibull 分布的特例。自动化的最佳拟合程序倾向于选择更灵活的分布（伽马和 Weibull 优于厄兰和指数），因为额外的灵活性使其能更贴合数据并获得更好的拟合汇总度量。但同样地，与数据贴合得紧密并不总能得到最恰当的输入模型。
3. 像 p 值这样的汇总统计量，终究只是一种汇总度量。它几乎不能说明拟合不佳出现在何处（在分布的中部、右尾还是左尾）。而人借助图形工具，能够看出拟合不佳出现在何处，并判断它对手头的应用是否重要。
我们的建议是：把自动化的分布选择当作提出候选分布的若干种方式之一。务必用图形方法检查自动选择的结果，并记住最终的选择权在你手中。
## 5 拟合非平稳泊松过程
一般而言，把非平稳泊松过程（nonstationary Poisson process，NSPP）拟合到到达数据是一个困难的问题，因为我们很少知道到达率函数 λ(t) 的恰当形式。一种途径是选择一个具有大量参数的非常灵活的模型，并用诸如最大似然之类的方法拟合它；关于这种途径的例子参见 Johnson、Lee 与 Wilson [1994]。第二种方法——也是这里我们考虑的方法——是把到达率近似为在某个基本时间区间（例如一小时、一天或一个月）内保持恒定，但在不同时间区间之间变化。这样，问题就变成选择基本时间区间以及估计每个区间内的到达率。
假设我们需要对某个时间段（比如 [0, T]）内的到达建模。我们所描述的方法最适用于这样一种情形：可以反复观测时间段 [0, T] 并对到达计数。例如，如果问题涉及对工作日（上午 8:00 至下午 6:00）期间电子邮件的到达建模，而我们相信到达率在半小时区间内近似恒定，那么我们需要能够统计若干天中各个半小时区间内的到达数。如果能够记录实际的到达时刻而非计数，那么实际到达时刻显然更好，因为事后可以把它们按我们希望的任何区间长度分组。不过，从此处往后我们将假定只有计数可用。
**表 8 NSPP 示例的星期一电子邮件到达数据**
| 时段 | 到达数 | 估计到达率（到达数/小时） | | |
|---|---|---|---|---|
| 第 1 天 | 第 2 天 | 第 3 天 | | |
| 8:00–8:30 | 12 | 14 | 10 | 24 |
| 8:30–9:00 | 23 | 26 | 32 | 54 |
| 9:00–9:30 | 27 | 19 | 32 | 52 |
| 9:30–10:00 | 20 | 13 | 12 | 30 |
把时间段 [0, T] 划分为 k 个长度为 \(\Delta t = T / k\) 的等长区间。例如，如果我们考虑从上午 8:00 到下午 6:00 的 10 小时工作日，并且允许速率每半小时变化一次，那么 \(T = 10 , k = 20\)，\(\Delta t = 1 / 2\)。在 n 个观测周期（例如 n 天）上，令 \(C_{i j}\) 为第 j 个观测周期的第 i 个时间区间内发生的到达数。在我们的例子中，\(C_{23}\) 就是第三个观测日从上午 8:30 到 9:00（第二个半小时区间）的到达数。
第 i 个时间段\(( i - 1 ) \Delta t < t \leq i \Delta t . ,\) 内的估计到达率，就是按时间区间长度缩放后的平均到达数：
$$
\widehat {\lambda} (t) = \frac {1}{n \Delta t} \sum_{j = 1} ^ {n} C_{i j}\tag{22}
$$
在各个时间区间的到达率估计出来之后，可以把速率看起来相同的相邻区间合并起来。
例如，考虑三个星期一工作日的前两个小时内的电子邮件到达计数，如表 8 所示。8:30–9:00 的估计到达率为
$$
\frac {1}{3 (1 / 2)} (23 + 26 + 32) = 54 \text{arrivals / hour}
$$
看到这些结果后，我们可能会考虑把区间 8:30–9:00 与区间 9:00–9:30 合并，因为两者的速率如此接近。还要注意，上一节描述的拟合优度检验可以分别应用于每个时间区间的数据，以检验泊松近似。
## 6 无数据时选择输入模型
遗憾的是，在实践中常常不得不在任何过程数据可用之前就建立仿真模型——也许是为了演示目的或做初步研究。在这种情况下，建模者必须在选择输入模型时足智多谋，并且必须仔细检查结果对所选模型的敏感性。
即使没有数据，也有若干途径可以获得关于某个过程的信息：
工程数据：产品或者过程常常带有制造商提供的性能额定值（例如，磁盘驱动器的平均失效时间为 10,000 小时；激光打印机每分钟可打印 8 页；刀具的切削速度为 1 厘米/秒；等等）。公司的规章可能规定时间或生产标准。这些数值通过确定一个中心值，为输入建模提供了起点。
专家意见：与对这个过程或类似过程有经验的人交谈。他们常常能够给出乐观的、悲观的和最可能的时间。他们也许还能说明这个过程是近乎恒定还是高度可变，并且可能能够指出变异性的来源。
物理或惯常限制：大多数真实过程在性能上存在物理限制；例如，计算机数据录入不可能比一个人打字的速度更快。由于公司政策的原因，一个过程可能允许耗费的时间存在上限。不要忽视那些缩小输入过程取值范围的显而易见的限制或界限。
过程的性质：第 2.2 节中对各种分布的描述，可以用来论证某个特定的选择，即使没有数据可用。
当数据不可得时，均匀分布、三角分布和 Beta 分布常被用作输入模型。均匀分布可能是糟糕的选择，因为在真实过程中，上限和下限很少与中心值一样可能。如果除了上限和下限之外还能给出一个最可能值，那么就可以使用三角分布。三角分布把大部分概率放在最可能值附近，而放在极端值附近的概率要小得多。如果使用 Beta 分布，那么一定要画出所选分布的密度函数；Beta 分布可能呈现不寻常的形状。
如果能够给出最小值、最大值以及一个或多个"断点"，就可以得到一种有用的改进。断点是一个中间值，连同小于或等于该值的概率。下面的例子说明如何使用断点。
### 例 20
对于一项生产计划仿真，需要各种产品的销售量。负责产品 XYZ-123 的销售人员说，销量不会少于 1000 件（因为已有合同），也不会多于 5000 件（因为那是该产品的整个市场）。根据她的经验，她认为销量超过 2000 件的机会是 90%，超过 3500 件的机会是 25%，而超过 4500 件的机会只有 1%。
表 9 汇总了这些信息。注意，超过某些销售目标的机会已被换算为小于或等于这些目标的累积概率。
**表 9 销售信息汇总**
| i | 区间（销售额） | 累积频率 \(c_i\) |
|---|---|---|
| 1 | \(1000 \le x \le 2000\) | 0.10 |
| 2 | \(2000 < x \le 3500\) | 0.75 |
| 3 | \(3500 < x \le 4500\) | 0.99 |
| 4 | \(4500 < x \le 5000\) | 1.00 |
当输入模型是在没有数据的情况下选定的时候，检验仿真结果对所选分布的敏感性就尤为重要。不仅要检验对分布中心的敏感性，还要检验对变异性或界限的敏感性。如果输出结果对输入模型极度敏感，这就提供了一个有说服力的理由，反对依据这些结果作出关键决策，而赞成着手开展数据收集工作。
## 7 多变量与时间序列输入模型
在第 1 节至第 4 节中，输入随机变量被认为与问题情境中的任何其他变量相互独立。当情况并非如此时，使用能够刻画这种依赖关系的输入模型就至关重要；否则，仿真可能极不准确。
### 例 21
某项供应链仿真包含工业机器人的提前期与年需求量。需求增加会导致提前期增加：这些机器人的最终组装必须按采购方的规格进行。因此，与其把提前期和需求当作相互独立的随机变量，不如建立一个多变量输入模型。
### 例 22
对某股票经纪人基于网络的交易站点的仿真，包含买入与卖出订单到达之间的时间间隔。投资者往往会对其他投资者的行为作出反应，因此这些买入与卖出订单是成批到达的。因此，与其把到达间隔时间当作相互独立的随机变量，不如建立一个时间序列模型。
我们把固定、有限个随机变量的多变量输入模型（例如例 21 中提前期与年需求这两个随机变量）与（概念上无限的）相关随机变量序列的时间序列输入模型（例如例 22 中相继的订单间隔时间）区分开来。在回顾两种依赖性度量——协方差与相关性——之后，我们将描述适用于这些例子的输入模型。
### 7.1 协方差与相关性
令 \(X_{1}\) 和 \(X_{2}\) 为两个随机变量，并令 \(\mu_{i} = \operatorname{E} ( X_{i} )\) 和 \(\sigma_{i} ^ {2} = \mathrm{V} ( X_{i} )\) 分别为 \(X_{i} ,\) 的均值和方差。协方差与相关性度量 \(X_{1}\) 与 \(X_{2} .\). 之间的线性依赖程度。换言之，协方差与相关性表明 \(X_{1}\) 与 \(X_{2}\) 之间的关系能被如下模型描述得多好：
$$
(X_{1} - \mu_{1}) = \beta (X_{2} - \mu_{2}) + \varepsilon
$$
其中 ε 是均值为 0 且独立于 \(X_{2}\) 的随机变量。如果事实上\(( X_{1} - \mu_{1} ) = \beta ( X_{2} - \mu_{2} )\)，那么这个模型就是完美的。另一方面，\(\operatorname{i f} X_{1}\) 与 \(X_{2}\) 统计独立，则 \(\beta = 0\)，该模型就没有价值。一般地，\(\beta\) 取正值表明 \(X_{1}\) 与 \(X_{2}\) 倾向于同时高于或同时低于各自的均值；\(\boldsymbol {\beta}\) 取负值表明它们往往位于各自均值的相反两侧。
\(X_{1}\) 与 \(X_{2}\) 之间的协方差定义为
$$
\operatorname{cov} \left(X_{1}, X_{2}\right) = \operatorname{E} \left[ \left(X_{1} - \mu_{1}\right) \left(X_{2} - \mu_{2}\right) \right] = \operatorname{E} \left(X_{1} X_{2}\right) - \mu_{1} \mu_{2}\tag{23}
$$
值 cov \(( X_{1} , X_{2} ) = 0\) 意味着在我们的依赖模型中 \(\beta = 0\)，而 \(\operatorname{c o v} ( X_{1} , X_{2} ) < 0 ( > 0 )\) 意味着 \(\beta < 0 ( > 0 )\)。协方差可以在 \(-\infty\) 与 \(\infty\) 之间取任意值。相关性把协方差标准化到介于 −1 与 1 之间：
$$
\rho = \operatorname{corr} (X_{1}, X_{2}) = \frac {\operatorname{cov} (X_{1} , X_{2})}{\sigma_{1} \sigma_{2}}\tag{24}
$$
同样，值 corr \(( X_{1} , X_{2} ) = 0\) 意味着在我们的模型中 \(\beta = 0\)，而 \(\operatorname{c o r r} ( X_{1} , X_{2} ) < 0 ( > 0 )\) 意味着 \(\beta < 0 ( > 0 )\)。\(\rho\) 越接近于 1 或 −1，\(X_{1}\) 与 \(X_{2}\) 之间的线性关系就越强。
现在，假设我们有一列随机变量 \(X_{1} , X_{2} , X_{3} , . . .\). 它们是同分布的（意味着它们都具有相同的均值和方差），但可能是相关的。我们把这样一列称为时间序列，并分别把 cov \(( X_{t} , X_{t + h} )\) 和 \(\operatorname{c o r r} ( X_{t} , X_{t + h} )\) 称为滞后 h 自协方差（lag-h autocovariance）和滞后 h 自相关（lag-h autocorrelation）。如果自协方差的取值只依赖于 h 而不依赖于 \(t ,\)，那么我们就说这个时间序列是协方差平稳的（covariance stationary）。对于协方差平稳的时间序列，我们使用如下简记
$$
\rho_{h} = \mathrm{corr} (X_{t}, X_{t + h})
$$
表示滞后 h 自相关。注意，自相关度量的是时间序列中相隔 \(h - 1\) 个其他变量的随机变量之间的依赖性。
### 7.2 多变量输入模型
如果 \(X_{1}\) 和 \(X_{2}\) 各自服从正态分布，那么它们之间的依赖性可以用参数为 \(\mu_{1} , \mu_{2} , \sigma_{1} ^ {2} , \sigma_{2} ^ {2}\) 和 \(\rho = \operatorname{c o r r} ( X_{1} , X_{2} )\)) 的二元正态分布来建模。\(\mu_{1} , \mu_{2} , \sigma_{1} ^ {2}\) 和 \(\sigma_{2} ^ {2}\) 的估计已在第 3.2 节中描述。为了估计 \(\rho ,\)，假设我们有 n 个独立同分布的对 \(( X_{11} , X_{21} ) , ( X_{12} , X_{22} ) , \dots , ( X_{1 n} , X_{2 n} )\)。那么样本协方差为
$$
\begin{array}{l} \widehat {\operatorname{cov}} (X_{1}, X_{2}) = \frac {1}{n - 1} \sum_{j = 1} ^ {n} \left(X_{1 j} - \bar {X} _ {1}\right) (X_{2 j} - \bar {X} _ {2}) \\ = \frac {1}{n - 1} \left(\sum_{j = 1} ^ {n} X_{1 j} X_{2 j} - n \bar {X} _ {1} \bar {X} _ {2}\right) \end{array}\tag{25}
$$
其中 \({\bar {X}} _ {1}\) 和 \({\bar {X}} _ {2}\) 是样本均值。相关性由下式估计
$$
\widehat {\rho} = \frac {\widehat {\mathrm{cov}} (X_{1} , X_{2})}{\widehat {\sigma} _ {1} \widehat {\sigma} _ {2}}\tag{26}
$$
其中 \(\widehat {\sigma} _ {1}\) 和 \(\widehat {\sigma} _ {2}\) 是样本方差。
### 例 23：例 21（续）
令 \(X_{1}\) 表示交货的平均提前期（以月计），\(X_{2}\) 表示工业机器人的年需求量。以下是过去十年关于需求与提前期的可得数据：
| 提前期 | 需求 |
|---|---|
| 6.5 | 103 |
| 4.3 | 83 |
| 6.9 | 116 |
| 6.0 | 97 |
| 6.9 | 112 |
| 6.9 | 104 |
| 5.8 | 106 |
| 7.3 | 109 |
| 4.5 | 92 |
| 6.3 | 96 |
常规计算给出 \({\bar {X}} _ {1} = 6 . 14 , {\widehat {\sigma}} _ {1} = 1 . 02 , {\bar {X}} _ {2} = 101 . 80 . ,\) 和 \(\widehat {\sigma} _ {2} = 9 . 93\)，分别作为 \(\mu_{1} , \sigma_{1} , \mu_{2}\) 和 \(\sigma_{2} .\) 的估计。为了估计相关性，我们需要
$$
\sum_{j = 1} ^ {10} X_{1 j} X_{2 j} = 6328. 5
$$
因此，\(\widehat {\mathrm{c o v}} = [ 6328 . 5 - ( 10 ) ( 6 . 14 ) ( 101 . 80 ) ] / ( 10 - 1 ) = 8 . 66\)，且
$$
\widehat {\rho} = \frac {8 . 66}{(1 . 02) (9 . 93)} = 0. 86
$$
显然，提前期与需求是强相关的。然而，在我们接受这个模型之前，还应当分别检验提前期和需求是否能被正态分布很好地表示。特别是，需求是一个离散取值的量，因此连续的正态分布充其量只是一种近似。
下面的简单算法可用于生成二元正态随机变量：
步骤 1. 生成 \(Z_{1}\) 和 \(Z_{2} . ,\)，即相互独立的标准正态随机变量。
步骤 2. 令 \(X_{1} = \mu_{1} + \sigma_{1} Z_{1}\)
步骤 3. 令 \(X_{2} = \mu_{2} + \sigma_{2} \left( \rho Z_{1} + \sqrt {1 - \rho^{2}} Z_{2} \right)\)
显然，二元正态分布并不适用于所有的多变量输入建模问题。它可以推广到 k 元正态分布，以建模两个以上随机变量之间的依赖性；但在许多情况下，任何形式的正态分布都是不恰当的。我们在第 7.4 节给出一种处理非正态分布的方法。关于其他模型的良好参考文献是 Johnson [1987] 以及 Nelson 与 Yamnitsky [1998]。
### 7.3 时间序列输入模型
如果 \(X_{1} , X_{2} , X_{3} , \dots\). 是一列同分布、但相关且协方差平稳的随机变量，那么有若干时间序列模型可用于表示这一过程。我们将描述两个具有如下特征的模型：其自相关具有形式
$$
\rho_{h} = \mathrm{corr} (X_{t}, X_{t + h}) = \rho^{h}
$$
其中 \(h = 1 , 2 , \ldots\)。注意 \(l a g .\)-h 自相关随着滞后增大而几何衰减，因此在时间上相距很远的观测值近乎独立。下面将要给出的一个模型中，每个 \(X_{t}\) 服从正态分布；在另一个模型中，每个 \(X_{t}\) 服从指数分布。更一般的时间序列输入模型在第 7.4 节以及 Nelson 与 Yamnitsky [1998] 中描述。
### AR(1) 模型
考虑时间序列模型
$$
X_{t} = \mu + \phi (X_{t - 1} - \mu) + \varepsilon_{t}\tag{27}
$$
其中 \(t = 2 , 3 , \ldots\)，\(\varepsilon_{2} , \varepsilon_{3} , \ldots .\). 独立同（正态）分布，均值为 0，方差为 \(\sigma_{\varepsilon} ^ {2}\)，且 \(- 1 < \phi < 1\)。如果初值 \(X_{1}\) 选择得当（见下文），那么 \(X_{1} , X_{2} , . . .\). 都服从均值为 \(\mu , ,\)、方差为 \(\sigma_{\varepsilon} ^ {2} / ( 1 - \phi^{2} )\) 的正态分布，且
$$
\rho_{h} = \phi^{h}
$$
其中 \(h = 1 , 2 , \ldots\)。这个时间序列模型称为 1 阶自回归模型（autoregressive order-1 model），简称 AR(1)。
参数 \(\phi\) 的估计可以由如下事实得到
$$
\phi = \rho^{1} = \mathrm{corr} (X_{t}, X_{t + 1})
$$
即滞后 1 自相关。因此，为了估计 \(\phi ,\)，我们首先用下式估计滞后 1 自协方差
$$
\widehat {\mathrm{cov}} (X_{t}, X_{t + 1}) = \frac {1}{n - 1} \sum_{t = 1} ^ {n - 1} (X_{t} - \bar {X}) (X_{t + 1} - \bar {X})
$$
$$
\doteq \frac {1}{n - 1} \left(\sum_{t = 1} ^ {n - 1} X_{t} X_{t + 1} - (n - 1) \bar {X} ^ {2}\right)\tag{28}
$$
并用通常的估计量 \({\widehat {\sigma}} ^ {2}\) 估计方差 \(\sigma^{2} = \operatorname{v a r} ( X )\)。然后
$$
\widehat {\phi} = \frac {\widehat {\mathrm{cov}} (X_{t} , X_{t + 1})}{\widehat {\sigma} ^ {2}}
$$
最后，分别用 \({\widehat {\mu}} = {\bar {X}}\) 和 \(\widehat {\sigma} _ {\varepsilon} ^ {2} = \widehat {\sigma} ^ {2} ( 1 - \widehat {\phi} ^ {2} )\)) 估计 \(\mu\) 和 \(\sigma_{\varepsilon} ^ {2}\)。
下面的算法在给定参数 \(\phi , \mu\) 和 \(\sigma_{\varepsilon} ^ {2} \mathrm{:}\) 的取值时，生成一个平稳的 \(\mathbf{A R} ( 1 )\)) 时间序列：
步骤 1. 从均值为 \(\mu\)、方差为 \(\sigma_{\varepsilon} ^ {2} / ( 1 - \phi^{2} )\) 的正态分布中生成 \(X_{1}\)。令 \(t = 2\)。
步骤 2. 从均值为 \(_ 0\)、方差为 \(\sigma_{\varepsilon} ^ {2}\) 的正态分布中生成 \(\varepsilon_{t}\)。
步骤 3. 令 \(X_{t} = \mu + \phi ( X_{t - 1} - \mu ) + \varepsilon_{t}\)。
步骤 4. 令 \(t = t + 1\)，转到步骤 2。
### EAR(1) 模型
考虑时间序列模型
$$
X_{t} = \left\{\begin{array}{l l} \phi X_{t - 1}, & \text{with probability} \phi \\ \phi X_{t - 1} + \varepsilon_{t}, & \text{with probability} 1 - \phi \end{array} \right.\tag{29}
$$
其中 \(t = 2 , 3 , \ldots\)，\(\varepsilon_{2} , \varepsilon_{3} , . . .\). 独立同（指数）分布，均值为 \(1 / \lambda\)，且 \(0 \leq \phi < 1\)。如果初值 \(X_{1}\) 选择得当（见下文），那么 \(X_{1} , X_{2} , \dots\) 都服从均值为 \(1 / \lambda\) 的指数分布，且
$$
\rho_{h} = \phi^{h}
$$
其中 \(h = 1 , 2 , \ldots\)。这个时间序列模型称为指数自回归 1 阶模型（exponential autoregressive order-1 model），简称 \(\mathrm{E A R} ( 1 )\)。该模型只能表示大于 0 的自相关。参数的估计按 \(\mathbf{A R} ( 1 )\) 的方式进行，即令 \({\widehat {\phi}} = {\widehat {\rho}} . ,\)（估计的滞后 1 自相关），并令 \(\widehat {\lambda} = 1 / \bar {X}\)。
下面的算法在给定参数 φ 和 λ 的取值时，生成一个平稳的 EAR(1) 时间序列：
步骤 1. 从均值为 \(1 / \lambda\) 的指数分布中生成 \(X_{1}\)。令 \(t = 2\)。
步骤 2. 从 [0, 1] 上的均匀分布中生成 U。如果 \(U \leq \phi .\)，则令
$$
X_{t} = \phi X_{t - 1}
$$
否则，从均值为 \(1 / \lambda\) 的指数分布中生成 \(\varepsilon_{t}\)，并令
$$
X_{t} = \phi X_{t - 1} + \varepsilon_{t}
$$
步骤 3. 令 \(t = t + 1\)，转到步骤 2。
### 例 24：例 22（续）
这位股票经纪人通常会有大量数据样本，但为了便于说明，假设已经记录了下列 20 个顾客买入与卖出订单之间的时间间隔（以秒计）：1.95、1.75、1.58、1.42、1.28、1.15、1.04、0.93、0.84、0.75、0.68、0.61、11.98、10.79、9.71、14.02、12.62、11.36、10.22、9.20。常规计算给出 \(\bar {X} = 5 . 2\) 和 \({\widehat {\sigma}} ^ {2} = 26 . 7\)。为了估计滞后 1 自相关，我们需要
$$
\sum_{j = 1} ^ {19} X_{t} X_{t + 1} = 924. 1
$$
于是，\({\widehat {\mathsf {c o v}}} = [ 924 . 1 - ( 20 - 1 ) ( 5 . 2 ) ^ {2} ] / ( 20 - 1 ) = 21 . 6 ,\)，且
$$
\widehat {\rho} = \frac {21 . 6}{26 . 7} = 0. 8
$$
因此，我们可以把到达间隔时间建模为 EAR(1) 过程，其中 \(\widehat {\lambda} = 1 / 5 . 2 = 0 . 192\) 和 \(\widehat {\phi} = 0 . 8 .\)，前提是各个间隔本身能被指数分布很好地建模。
### 7.4 正态到任意分布变换
二元正态分布以及 AR(1) 和 EAR(1) 时间序列模型都是有用的输入模型，它们既容易拟合也容易仿真。然而，其边缘分布要么是正态的，要么是指数的，这对许多应用来说当然不是最好的选择。幸运的是，我们可以从二元正态或 AR(1) 模型出发，把它变换成具有我们所想要的任意边缘分布（包括指数分布）。
假设我们要仿真一个 cdf 为 \(F ( x )\) 的随机变量 X。令 \(Z\) 为标准正态随机变量（均值为 0，方差为 1），\(\Phi ( z )\) 为其 cdf。那么可以证明
$$
R = \Phi (Z)
$$
是一个 \(U ( 0 ,\)1) 随机变量。如果我们有一个 \(U ( 0 , 1 )\) 随机变量，就可以通过逆 cdf 变换得到 X：
$$
X = F^{- 1} [ R ] = F^{- 1} [ \Phi (Z) ]
$$
我们称之为正态到任意分布变换（normal to anything transformation），简称 NORTA。
当然，如果我们想要的只是 \(X ,\)，那就没有理由如此周折；我们可以直接生成 R。但假设我们想要一个二元随机向量 \(( X_{1} , X_{2} )\)，使得 \(X_{1}\) 与 \(X_{2}\) 相关，但它们的分布不是正态的。那么我们可以从二元正态随机向量 \(( Z_{1} , Z_{2} )\) 出发，应用 NORTA 变换得到
$$
X_{1} = F_{1} ^ {- 1} [ \Phi (Z_{1}) ] \text{and} X_{2} = F_{2} ^ {- 1} [ \Phi (Z_{2}) ]
$$
甚至并不要求 \(F_{1}\) 与 \(F_{2}\) 来自同一个分布族；例如，\(F_{1}\) 可以是某个指数分布，而 \(F_{2}\) 是某个 Beta 分布。
同样的思想也适用于时间序列。如果 \(Z_{t}\) 由边缘分布为 \(N ( 0 , 1 )\) 的 \(\mathbf{A R} ( 1 )\) 生成，那么
$$
X_{t} = F^{- 1} [ \Phi (Z_{t}) ]
$$
就是一个边缘分布为 \(F ( x )\) 的时间序列模型。为了确保 \(Z_{t}\) 服从 \(N ( 0 , 1 )\)，我们在 AR(1) 模型中令 \(\mu = 0\) 和 \(\sigma _ {\mathrm {{s}}} ^ {2} = 1 - \phi^{2}\)。
尽管 NORTA 方法非常通用，但要实施它必须解决两个技术问题：
1. NORTA 方法需要能够计算标准正态 cdf\(\Phi ( z )\)，以及相关分布的逆 cdf\({\cal F} ^ {- 1} ( R )\)。\(\Phi ( z )\) 没有闭式表达式，许多分布的 \({\cal F} ^ {- 1} ( R )\) 也没有闭式表达式。因此需要数值近似。幸运的是，这些函数已内置在许多符号计算和电子表格程序中，下面我们就给出一个例子。此外，Bratley、Fox 与 Schrage [1987] 收录了针对许多分布的算法。
2. 标准正态随机变量 \(( Z_{1} , Z_{2} )\) 之间的相关性在通过 NORTA 变换时会被扭曲。更具体地说，如果 \(( Z_{1} , Z_{2} )\)) 的相关系数为 \(\rho ,\)，那么一般地 \(X_{1} = F_{1} ^ {- 1} [ \Phi ( Z_{1} ) ]\) 与 \(X_{2} = F_{2} ^ {- 1} [ \Phi ( Z_{2} ) ]\)] 将具有相关系数 \(\rho_{X} \neq \rho\)。这一差异常常很小，但并非总是如此。
第二个问题更为关键，因为在输入建模问题中，我们希望指定二元相关或滞后 1 相关。因此，我们需要找到能给出我们所想要的输入相关 \(\rho_{X}\) 的二元正态相关系数 \(\rho\)（回顾一下，我们通过滞后 1 相关\(\rho_{X} = \operatorname{c o r r} ( X_{t} , X_{t + 1} ) )\) 来指定时间序列模型）。关于这一问题已有大量研究，包括 Cario 与 Nelson [1996, 1998] 以及 Biller 与 Nelson [2003]。幸运的是，已经证明 \(\rho_{X}\) 是
关于 \(\rho ,\) 的非降函数，且 \(\rho\) 与 \(\rho_{X}\) 总是同号。因此，我们可以基于如下算法做一个相对简单的搜索：
步骤 1. 先令 \(\rho = \rho_{X}\)。
步骤 2. 生成大量的相关系数为 \(\rho , ,\) 的二元正态对 \(( Z_{1} , Z_{2} )\))，并用 NORTA 变换把它们变换为 \(( X_{1} , X_{2} ) \mathbf{s}\)。
步骤 3. 用式 (25) 计算 \(( X_{1} , X_{2} )\) 之间的样本相关，记为 \(\widehat {\rho} _ {T}\)。如果 \(\widehat {\rho} _ {T} > \rho_{X}\)，则减小 \(\rho\) 并转到步骤 2；如果 \(\widehat {\rho_{T}} < \rho_{X}\)，则增大 \(\rho\) 并 \({\bf g 0}\)转到步骤 2。如果 \(\widehat {\rho} _ {T} \approx \rho_{X}\)，则停止。
### 例 25
假设我们需要 \(X_{1}\) 服从均值为 1 的指数分布，\(X_{2}\) 服从 \(\beta_{1} = 1 , \beta_{2} = 1 / 2\) 的 Beta 分布，并且二者具有相关系数 \(\rho_{X} = 0 . 45\)。图 8 给出一个用 Maple 编写的程序，它将估计所需的 \(\rho .\) 值。在该程序中，n 是用于估计相关系数的样本对数。取 n = 1000 运行该程序，得到 \(\rho = 0 . 52\)。
```txt
NORTARho := proc(rhoX, n)
local Z1, Z2, Ztemp, X1, X2, R1, R2, rho, rhoT, lower, upper;
randomize(123456);
Z1 := [random[normald[0,1]](n)]:
ZTemp := [random[normald[0,1]](n)]:
Z2 := [0]:
# set up bisection search
rho := rhoX:
if (rhoX < 0) then
lower := -1:
upper := 0:
else
lower := 0:
upper := 1:
fi:
Z2 := rho*Z1 + sqrt(1-rho^2)*ZTemp:
R1 := statevalf[cdf,normald[0,1]](Z1):
R2 := statevalf[cdf,normald[0,1]](Z2):
X1 := statevalf/icdf,exponential[1,0]](R1):
X2 := statevalf/icdf,beta[1,2]](R2):
rhoT := describe(linearcorrelation)(X1, X2);
# do bisection search until 5% relative error
while abs(rhoT - rhoX)/abs(rhoX) > 0.05 do
if (rhoT > rhoX) then
upper := rho:
else
lower := rho:
fi:
rho := evalf((lower + upper)/2):
Z2 := rho*Z1 + sqrt(1-rho^2)*ZTemp:
R1 := statevalf[cdf,normald[0,1]](Z1):
R2 := statevalf[cdf,normald[0,1]](Z2):
X1 := statevalf/icdf,exponential[1,0]](R1):
X2 := statevalf/icdf,beta[1,2]](R2):
rhoT := describe(linearcorrelation)(X1, X2);
end do;
RETURN(rho);
end;
```
**图 8 估计 NORTA 方法所需的二元正态相关系数的 Maple 程序**
## 8 小结
在离散事件仿真项目中，输入数据的收集与分析需要投入大量的时间与资源。然而，无论仿真模型多么有效、多么精巧，不可靠的输入都会导致输出在后续解释时可能产生错误的建议。
本章讨论了建立输入数据模型的四个步骤：收集原始数据、识别底层的统计分布、估计参数以及检验拟合优度。
针对便利数据收集工作给出了一些建议。然而，诸如完成习题 1 至习题 5 中任何一道所获得的经验，将会增强人们对数据收集中可能出现的困难以及规划之必要性的认识。
数据一旦收集完毕，就应当假设一个统计模型。此时如果有足够的数据，绘制直方图非常有用。通常可以基于底层过程以及直方图的形状选择一个分布，作进一步研究。
研究继续进行，对所假设的分布估计参数。本章给出了仿真中常用分布的建议估计量。在许多情况下，这些估计量是样本均值与样本方差的函数。
该过程的最后一步是检验分布假设。\(q {-} q\) 图是评估拟合的有用图形方法。卡方检验、Kolmogorov–Smirnov 检验和 Anderson-Darling 拟合优度检验可应用于许多分布假设。当某个分布假设被拒绝时，就尝试另一个分布。当所有办法都失败时，可以在模型中使用经验分布。
遗憾的是，在某些情况下，必须在没有时间或资源去收集作为输入模型依据的数据时开展仿真研究。当这种情况发生时，分析人员必须利用任何可得的信息——例如制造商规格和专家意见——来构造输入模型。当输入模型是在没有数据支持的情况下导出时，检查结果对所选模型的敏感性就尤为重要。
许多（但并非全部）输入过程都可以表示为独立同分布随机变量的序列。当输入应当表现出依赖性时，就需要多变量输入模型。二元正态分布（更一般地，多元正态分布）常被用来表示有限个相关的随机变量。时间序列模型则适用于表示（概念上无限的）相关输入序列。NORTA 变换便于建立边缘分布非正态的多变量输入模型。
## 参考文献
BILLER, B., AND B. L. NELSON [2003], "Modeling and Generating Multivariate Time Series with Arbitrary Marginals Using an Autoregressive Technique," *ACM Transactions on Modeling and Computer Simulation*, Vol. 13, pp. 211–237.
BRATLEY, P., B. L. FOX, AND L. E. SCHRAGE [1987], *A Guide to Simulation*, 2d ed., Springer-Verlag, New York.
CARIO, M. C., AND B. L. NELSON [1996], "Autoregressive to Anything: Time-Series Input Processes for Simulation," *Operations Research Letters*, Vol. 19, pp. 51–58.
CARIO, M. C., AND B. L. NELSON [1998], "Numerical Methods for Fitting and Simulating Autoregressive-to-Anything Processes," *INFORMS Journal on Computing*, Vol. 10, pp. 72–81.
CHOI, S. C., AND R. WETTE [1969], "Maximum Likelihood Estimation of the Parameters of the Gamma Distribution and Their Bias," *Technometrics*, Vol. 11, No. 4, pp. 683–890.
CHAMBERS, J. M., W. S. CLEVELAND, AND P. A. TUKEY [1983], *Graphical Methods for Data Analysis*, CRC Press, Boca Raton, FL.
CONNOVER, W. J. [1998], *Practical Nonparametric Statistics*, 3d ed., Wiley, New York.
DURBIN, J. [1975], "Kolmogorov–Smirnov Tests When Parameters Are Estimated with Applications to Tests of Exponentiality and Tests on Spacings," *Biometrika*, Vol. 65, pp. 5–22.
FISHMAN, G. S. [1973], *Concepts and Methods in Discrete Event Digital Simulation*, Wiley, New York.
GUMBEL, E. J. [1943], "On the Reliability of the Classical Chi-squared Test," *Annals of Mathematical Statistics*, Vol. 14, pp. 253ff.
HINES, W. W., D. C. MONTGOMERY, D. M. GOLDSMAN, AND C. M. BORROR [2002], *Probability and Statistics in Engineering and Management Science*, 4th ed., Wiley, New York.
JOHNSON, M. A., S. LEE, AND J. R. WILSON [1994], "NPPMLE and NPPSIM: Software for Estimating and Simulating Nonhomogeneous Poisson Processes Having Cyclic Behavior," *Operations Research Letters*, Vol. 15, pp. 273–282.
JOHNSON, M. E. [1987], *Multivariate Statistical Simulation*, Wiley, New York.
LAW, A. M. [2007], *Simulation Modeling and Analysis*, 4th ed., McGraw-Hill, New York.
LILLIEFORS, H. W. [1967], "On the Kolmogorov–Smirnov Test for Normality with Mean and Variance Unknown," *Journal of the American Statistical Association*, Vol. 62, pp. 339–402.
LILLIEFORS, H. W. [1969], "On the Kolmogorov–Smirnov Test for the Exponential Distribution with Mean Unknown," *Journal of the American Statistical Association*, Vol. 64, pp. 387–389.
MANN, H. B., AND A. WALD [1942], "On the Choice of the Number of Intervals in the Application of the Chi-squared Test," *Annals of Mathematical Statistics*, Vol. 18, p. 50ff.
NELSON, B. L., AND M. YAMNITSKY [1998], "Input Modeling Tools for Complex Problems," in *Proceedings of the 1998 Winter Simulation Conference*, D. Medeiros, E. Watson, J. Carson, and M. Manivannan, eds., Washington, DC, Dec. 13–16, pp. 105–112.
STUART, A., J. K. ORD, AND E. ARNOLD [1998], *Kendall's Advanced Theory of Statistics*, 6th ed., Vol. 2, Oxford University Press, Oxford.
## 习题
1. 去一家小商店，记录到达间隔时间与服务时间的分布。如果有若干名工作人员，各人的服务时间分布彼此相比如何？是否需要为每一类器具构造服务时间分布？（务必取得管理层许可后再开展此项研究。）
2. 去一家自助餐厅，收集到达间隔时间与服务时间分布的数据。到达间隔时间分布很可能对一天三餐各不相同，并且在一餐之内也可能变化——也就是说，上午 11:00 至中午 12:00 的到达间隔时间分布可能与中午 12:00 至下午 1:00 的不同。把服务时间定义为：从顾客到达可以作出第一项选择的位置之时起，到其离开餐厅队列之时为止。（对这一定义作任何合理的修改都是可以接受的。）服务时间分布很可能对每一餐都不同。就这两种分布而言，能否把一天中的某些时段或一周中的某些日子分组，以适应数据的同质性？（务必取得管理层许可后再开展此项研究。）
3. 去一个主要交通交叉口，记录来自各个方向的到达间隔时间分布。有些到达者要直行，有些要左转，有些要右转。到达间隔时间分布在一天之中变化，也随一周中的星期几而变化。时不时还会发生事故。为你的数据建立一个输入模型。
4. 去一家杂货店，构造收银台处的到达分布与服务分布。这些分布可能随一天中的时段和一周中的星期几而变化。同时记录所有时间可用的服务通道数量。（务必取得管理层许可后再开展此项研究。）
5. 去一家自助洗衣店，收集数据以建立"一位顾客使用的洗衣机与烘干机数量"（很可能相关）的输入模型。同时收集顾客到达率的数据（很可能不是平稳的）。（务必取得管理层许可后再开展此项研究。）
6. 在同一张图上绘制四个理论正态密度函数，每个分布的均值都为零，但标准差分别为 1/4、1/2、1 和 2。
7. 在同一张图上绘制厄兰分布的 pdf，其中 \(\theta = 1 / 2\) 且 \(k = 1 , 2 , 4\), 和 8。
8. 在同一张图上绘制厄兰分布的 pdf，其中 θ = 2 且 k = 1、2、4 和 8。
9. 绘制参数 α 取下列各值时泊松分布的 pmf：
$$
\alpha = 1 / 2
$$
$$
\alpha = 1
$$
$$
\alpha = 2
$$
$$
\alpha = 4
$$
10. 在同一张图上绘制参数 λ 分别等于 0.6 和 1.2 时得到的两条指数分布 pdf。
11. 在同一张图上绘制 \(\nu = 0 , \alpha = 1 / 2 ,\) 且 \(\beta = 1 , 2 ,\) 和 4 时得到的三条 Weibull 分布 pdf。
12. 下列数据是从一个伽马分布随机生成的：
| 1.691 | 1.437 | 8.221 | 5.976 |
|---|---|---|---|
| 1.116 | 4.435 | 2.345 | 1.782 |
| 3.810 | 4.589 | 5.313 | 10.90 |
| 2.649 | 2.432 | 1.581 | 2.432 |
| 1.843 | 2.466 | 2.833 | 2.361 |
计算最大似然估计量 \(\widehat {\beta}\) 和 \({\widehat {\theta}} .\)。
13. 下列数据是从一个 \(\nu = 0 \mathrm{:}\) 的 Weibull 分布随机生成的：
| 7.936 | 5.224 | 3.937 | 6.513 |
|---|---|---|---|
| 4.599 | 7.563 | 7.172 | 5.132 |
| 5.259 | 2.759 | 4.278 | 2.696 |
| 6.212 | 2.407 | 1.857 | 5.002 |
| 4.612 | 2.003 | 6.908 | 3.326 |
计算最大似然估计量 α 和 \({\widehat {\beta}} .\)。（本题需要一台可编程计算器、一台计算机，或者极大的耐心。）
14. 佐治亚州亚特兰大市与雅典市之间的公路沿线 100 公里内事故高发。公共安全官员称，沿线事故的发生是随机（均匀）分布的，但新闻媒体则不这么认为。佐治亚州公共安全部公布了 9 月份的记录。这些记录给出了 30 起涉及伤亡的事故发生的地点如下（数据点表示距亚特兰大市界的距离）：
| 88.3 | 40.7 | 36.3 | 27.3 | 36.8 |
|---|---|---|---|---|
| 91.7 | 67.3 | 7.0 | 45.2 | 23.3 |
| 98.8 | 90.1 | 17.2 | 23.7 | 97.4 |
| 32.4 | 87.8 | 69.8 | 62.6 | 99.7 |
| 20.6 | 73.1 | 21.6 | 6.0 | 45.3 |
| 76.6 | 73.2 | 27.3 | 87.6 | 87.2 |
用 Kolmogorov–Smirnov 检验判断 9 月份事故地点的分布是否服从均匀分布。
15. 证明例 19 的 Kolmogorov–Smirnov 检验统计量为 \(D = 0 . 1054\)。
16. 某联邦机构正在研究一座地下煤矿每月工伤事故数量的记录。过去 100 个月的数值如下：
| 每月工伤数 | 出现频数 |
|---|---|
| 0 | 35 |
| 1 | 40 |
| 2 | 13 |
| 3 | 6 |
| 4 | 4 |
| 5 | 1 |
| 6 | 1 |
(a) 对这些数据应用卡方检验，检验"底层分布为泊松分布"的假设。取显著性水平 \(\alpha = 0 . 05\)。
(b) 对这些数据应用卡方检验，检验"该分布是均值为 1.0 的泊松分布"的假设。同样取 \(\alpha = 0 . 05 .\)。
(c) (a) 与 (b) 两部分有何区别，各自可能在什么情形下出现？
17. 对 50 名不同员工计算和记录本周工作小时数所需的时间进行了测量，结果如下（以分钟计）：
| 员工 | 时间（分钟） | 员工 | 时间（分钟） |
|---|---|---|---|
| 1 | 1.88 | 26 | 0.04 |
| 2 | 0.54 | 27 | 1.49 |
| 3 | 1.90 | 28 | 0.66 |
| 4 | 0.15 | 29 | 2.03 |
| 5 | 0.02 | 30 | 1.00 |
| 6 | 2.81 | 31 | 0.39 |
| 7 | 1.50 | 32 | 0.34 |
| 8 | 0.53 | 33 | 0.01 |
| 9 | 2.62 | 34 | 0.10 |
| 10 | 2.67 | 35 | 1.10 |
| 11 | 3.53 | 36 | 0.24 |
| 12 | 0.53 | 37 | 0.26 |
| 13 | 1.80 | 38 | 0.45 |
| 14 | 0.79 | 39 | 0.17 |
| 15 | 0.21 | 40 | 4.29 |
| 16 | 0.80 | 41 | 0.80 |
| 17 | 0.26 | 42 | 5.50 |
| 18 | 0.63 | 43 | 4.91 |
| 19 | 0.36 | 44 | 0.35 |
| 20 | 2.03 | 45 | 0.36 |
| 21 | 1.42 | 46 | 0.90 |
| 22 | 1.28 | 47 | 1.03 |
| 23 | 0.82 | 48 | 1.73 |
| 24 | 2.16 | 49 | 0.38 |
| 25 | 0.05 | 50 | 0.48 |
用卡方检验（如例 18 那样）检验"这些服务时间服从指数分布"的假设。取类别区间数 k = 6。取显著性水平 α = 0.05。
18. Studentwiser Beer Company 正在设法查明其玻璃瓶的抗破强度分布。随机抽取 50 个瓶子并测试其抗破强度，结果如下（单位：磅/平方英寸）：
| 218.95 | 232.75 | 212.80 | 231.10 | 215.95 |
|---|---|---|---|---|
| 237.55 | 235.45 | 228.25 | 218.65 | 212.80 |
| 230.35 | 228.55 | 216.10 | 229.75 | 229.00 |
| 199.75 | 225.10 | 208.15 | 213.85 | 205.45 |
| 219.40 | 208.15 | 198.40 | 238.60 | 219.55 |
| 243.10 | 198.85 | 224.95 | 212.20 | 222.90 |
| 218.80 | 203.35 | 223.45 | 213.40 | 206.05 |
| 229.30 | 239.20 | 201.25 | 216.85 | 207.25 |
| 204.85 | 219.85 | 226.15 | 230.35 | 211.45 |
| 227.95 | 229.30 | 225.25 | 201.25 | 216.10 |
使用输入建模软件，应用软件中所有可用的正态性检验。如果有卡方检验可用，至少用两种不同的区间数选择来应用它。所有检验是否得出相同的结论？
19. Crosstowner 是一条从亚特兰大东北部斜穿到西南部的公交线路。完成该线路所需的时间由公交司机记录。该线路周一至周五运营。最近五十次上午 8:00 班次的运行时间如下（以分钟计）：
| 92.3 | 92.8 | 106.8 | 108.9 | 106.6 |
|---|---|---|---|---|
| 115.2 | 94.8 | 106.4 | 110.0 | 90.9 |
| 104.6 | 72.0 | 86.0 | 102.4 | 99.8 |
| 87.5 | 111.4 | 105.9 | 90.7 | 99.2 |
| 97.8 | 88.3 | 97.5 | 97.4 | 93.7 |
| 99.7 | 122.7 | 100.2 | 106.5 | 105.5 |
| 80.7 | 107.9 | 103.2 | 116.4 | 101.7 |
| 84.8 | 101.9 | 99.1 | 102.2 | 102.5 |
| 111.7 | 101.5 | 95.1 | 92.8 | 88.5 |
| 74.4 | 98.9 | 111.9 | 96.5 | 95.9 |
这些运行时间服从什么分布？建立并检验一个合适的模型。
20. 在某通信中心以电子方式对一条报文的传输时间（以毫秒计）进行抽样。样本中最后 50 个取值如下：
| 7.936 | 4.612 | 2.407 | 4.278 | 5.132 |
|---|---|---|---|---|
| 4.599 | 5.224 | 2.003 | 1.857 | 2.696 |
| 5.259 | 7.563 | 3.937 | 6.908 | 5.002 |
| 6.212 | 2.759 | 7.172 | 6.513 | 3.326 |
| 8.761 | 4.502 | 6.188 | 2.566 | 5.515 |
| 3.785 | 3.742 | 4.682 | 4.346 | 5.359 |
| 3.535 | 5.061 | 4.629 | 5.298 | 6.492 |
| 3.502 | 4.266 | 3.129 | 1.298 | 3.454 |
| 5.289 | 6.805 | 3.827 | 3.912 | 2.969 |
| 4.646 | 5.963 | 3.829 | 4.404 | 4.924 |
这些传输时间服从什么分布？建立并检验一个合适的模型。
21. 在 Gotwatts Flash and Flicker Company 精确记录了电力服务接通申请之间的时间间隔（以分钟计），最近 50 次申请的结果如下：
| 0.661 | 4.910 | 8.989 | 12.801 | 20.249 |
|---|---|---|---|---|
| 5.124 | 15.033 | 58.091 | 1.543 | 3.624 |
| 13.509 | 5.745 | 0.651 | 0.965 | 62.146 |
| 15.512 | 2.758 | 17.602 | 6.675 | 11.209 |
| 2.731 | 6.892 | 16.713 | 5.692 | 6.636 |
| 2.420 | 2.984 | 10.613 | 3.827 | 10.244 |
| 6.255 | 27.969 | 12.107 | 4.636 | 7.093 |
| 6.892 | 13.243 | 12.711 | 3.411 | 7.897 |
| 12.413 | 2.169 | 0.921 | 1.900 | 0.315 |
| 4.370 | 0.377 | 9.063 | 1.875 | 0.790 |
这些服务申请之间的时间间隔服从什么分布？建立并检验一个合适的模型。
22. Earth Moving Tractor Company 保存了 D-3 拉铲挖掘机变速箱大修套件的日需求记录，结果如下：
| 0 | 2 | 0 | 0 | 0 |
|---|---|---|---|---|
| 1 | 0 | 1 | 1 | 1 |
| 0 | 1 | 0 | 0 | 0 |
| 2 | 0 | 1 | 0 | 1 |
| 0 | 1 | 0 | 0 | 2 |
| 1 | 0 | 1 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 |
| 1 | 0 | 1 | 0 | 1 |
| 0 | 0 | 3 | 0 | 1 |
| 1 | 0 | 0 | 0 | 0 |
这些日需求服从什么分布？建立并检验一个合适的模型。
23. 将对一个执行两道工序（先铣削、再刨削）的作业车间开展仿真。本来可以分别收集两道工序加工时间的数据，然后从各自的分布中生成随机发生值。然而，车间经理说这两个时间可能是相关的；大的铣削作业需要大量刨削。对随后的 25 个订单收集了数据，结果如下（以分钟计）：
| 订单 | 铣削时间（分钟） | 刨削时间（分钟） | 订单 | 铣削时间（分钟） | 刨削时间（分钟） |
|---|---|---|---|---|---|
| 1 | 12.3 | 10.6 | 14 | 24.6 | 16.6 |
| 2 | 20.4 | 13.9 | 15 | 28.5 | 21.2 |
| 3 | 18.9 | 14.1 | 16 | 11.3 | 9.9 |
| 4 | 16.5 | 10.1 | 17 | 13.3 | 10.7 |
| 5 | 8.3 | 8.4 | 18 | 21.0 | 14.0 |
| 6 | 6.5 | 8.1 | 19 | 19.5 | 13.0 |
| 7 | 25.2 | 16.9 | 20 | 15.0 | 11.5 |
| 8 | 17.7 | 13.7 | 21 | 12.6 | 9.9 |
| 9 | 10.6 | 10.2 | 22 | 14.3 | 13.2 |
| 10 | 13.7 | 12.1 | 23 | 17.0 | 12.5 |
| 11 | 26.2 | 16.0 | 24 | 21.2 | 14.2 |
| 12 | 30.4 | 18.9 | 25 | 28.4 | 19.1 |
| 13 | 9.9 | 7.7 | | | |
(a) 以铣削时间为横轴、刨削时间为纵轴作图。这些数据看起来是否相关？
(b) 计算铣削时间与刨削时间之间的样本相关系数。
(c) 对这些数据拟合一个二元正态分布。
24. 编写一个计算机程序，计算 Weibull 分布的最大似然估计量 \(( {\widehat {\alpha}} , {\widehat {\beta}} )\)。程序的输入应当包括样本量 \(n ;\)、观测值 \(x_{1} , x_{2} , \ldots , x_{n} ;\)、停止准则 （当 $| f ( \widehat {\beta} _ {j} ) | \leq \epsilon ) $ 时停止）；以及打印选项 OPT（通常设为 0）。输出为估计值 α 和 \({\widehat {\beta}} . {\mathrm{I f}} {\mathrm{O P T}} = 1\)，将打印附加输出，如表 4 那样显示收敛过程。使程序尽可能"用户友好"。
25. 考察一个你可以访问的计算机软件库或仿真支持环境。获取其中对求解习题 7 至习题 24 有用的数据分析软件的文档。用该软件辅助求解选定的题目。
26. 某家小旅馆连续 20 晚的住客人数被观测为 20、14、21、19、14、18、21、25、27、26、22、18、13、18、18、18、25、23、20、21。对这些数据分别拟合 AR(1) 与 EAR(1) 模型。通过查看数据的直方图判断哪个模型拟合更好。
27. 下列数据表示在一家银行办理业务所需的时间，以分钟计：0.740、1.28、1.46、2.36、0.354、0.750、0.912、4.44、0.114、3.08、3.24、1.10、1.59、1.47、1.17、1.27、9.12、11.5、2.42、1.77。为这些数据建立一个输入模型。
28. 两类作业（A 与 B）随订单到达而被释放到某作业车间的输入缓冲区，订单的到达是不确定的。上周的生产数据如下：
| 天 | 作业数 | A 的数量 |
|---|---|---|
| 1 | 83 | 53 |
| 2 | 93 | 62 |
| 3 | 112 | 66 |
| 4 | 65 | 41 |
| 5 | 78 | 55 |
建立每天各类作业新到达数的输入模型。
29. 下列是某台机器上加工时间（以分钟计）的可得数据：0.64、0.59、1.1、3.3、0.54、0.04、0.45、0.25、4.4、2.7、2.4、1.1、3.6、0.61、0.20、1.0、0.27、1.7、0.04、0.34。建立加工时间的输入模型。
30. 对于此处给出的到达计数数据，估计 NSPP 的到达率函数 λ(t)，以每小时到达数计：
| 时段 | 到达数 | | | |
|---|---|---|---|---|
| 第 1 天 | 第 2 天 | 第 3 天 | 第 4 天 | |
| 8:00–10:00 | 22 | 24 | 20 | 28 |
| 10:00–12:00 | 23 | 26 | 32 | 30 |
| 12:00–2:00 | 40 | 33 | 32 | 38 |
31. 利用网络，调研本章提到的某些输入建模软件包（可以从 www.informs-sim.org 入手）。它们有哪些功能？包含哪些分布？
32. 在外科收治流程中，有护工负责护送患者前往术前准备区。按政策规定，患者不得自行前往。一位管理工程师快步走了一遍该路线，发现走到准备区需要 2 分钟。护工说，护送一位患者到该区域有时需要长达 7 分钟，但通常需要 5 分钟。根据所给信息，构造护送患者所需时间的输入模型。
33. 在一道插接/密封/检验作业中，电子元件被插装到电路板上，外壳被组装并密封，然后对成品进行测试。完成这一切的时间预计对 A12117c 型组件为 2 分钟，对 B33433x 型组件为 3 分钟。围绕这些值上下 20% 的平均偏差被视为标准的变异性容许量。根据所给信息，构造该插接/密封/检验作业的输入模型。
34. 下列数据是在五个星期一按一天中的小时统计的呼叫中心来话计数。该过程看起来是否是非平稳的？如果是，则为上午 8:00 至下午 4:00 之间呼叫的到达建立一个 NSPP 输入模型。
| 8-9A.M. | 9-10A.M. | 10-11A.M. | 11-12P.M. | 12-1P.M. | 1-2P.M. | 2-3P.M. | 3-4P.M. |
|---|---|---|---|---|---|---|---|
| 27 | 37 | 58 | 68 | 65 | 33 | 37 | 21 |
| 27 | 35 | 67 | 93 | 74 | 42 | 39 | 21 |
| 25 | 39 | 58 | 75 | 82 | 34 | 46 | 19 |
| 18 | 42 | 48 | 75 | 62 | 51 | 49 | 14 |
| 29 | 42 | 65 | 88 | 70 | 43 | 45 | 22 |