系统建模与仿真课件

随机变量的生成

进行随机仿真，必须要随机抽样，即产生服从一定分布的随机变量。若给出随机变量的分布函数F(x),，则可用各种方法生成服从该分布的随机变量。

1、逆变法

2、函数变换法

一、原理对任一严格单调增的分布函数F(x)，其逆函数X=F-1(U)的分布函数为若F(x)是随机变量X的严格单调的分布函数，其逆函数U=F(X)也是随机变量，0≤F(X)≤1，故0≤U≤1，于是对任一0≤u≤1有

逆变法

通过求随机变量的分布函数F(x)的逆函数，得到分布为F(x)的随机数。由于通过逆函数得到随机数，故称为逆变法。逆变法的步骤是：

①生成随机数；

②用逆函数产生随机变量。二、均匀分布

在区间[a,b]上均匀分布U(a,b)，其概率密度函数f(x)及分布函数F(x)分别是：

则其逆函数为用随机数发生器xi＝75xi-1mod(231-1)及Excel的随机数发生器产生随机数，据此试生成区间[3,5]服从均匀分布的随机变量，并比较两种结果。解：①随机数发生器xi＝75xi-1mod(231-1)，ui=xi/(231-1)和Excel的随机数发生器产生随机数ui；②由公式生成[3,5]的随机变量xi=3+(5-3)*ui，结果及随机变量与随机数的关系图如下。三、指数分布指数分布的概率密度函数及分布函数为指数分布的概率密度函数及分布函数为令u＝F-1(x)（u

为[0,1]均匀分布的随机数），即解此方程，有

生成均值为3的指数分布随机变量。解：利用Excel的随机数发生器产生[0,1]区间的随机数u，由逆变法公式

得指数分布随机变量，计算结果如下四、离散分布X是离散随机变量，取值为x1,x2,…,xn，其概率分布为分布函数为xF(x)x1x2x3

0

x4xnp1p2p3p4pn1离散分布随机变量的分布函数

计算离散分布随机变量的步骤是：

①产生U(0,1)随机数；

②若u≤p1，则x=x1；否则当

时，x=xi，x是分布函数为F(x)的随机数。用逆变法求离散分布的随机变量。解：先用随机数发生器xi＝75xi-1mod(231-1)产生随机数；比较所得随机数u与分布函数，若u≤p1，x=x1；否则当

时，x=xi，结果如下：某产品的单位成本随市场随机波动，历史数据统计如下：

产品单位成本分布概率试用逆变法产生随机变量。序号

单位成本（元）

概率1430.102440.203450.404460.205470.10逆变法须先求出分布函数的逆函数F-1，而有些分布函数的解析逆函数表达式难于求得，虽可用数值计算或幂级数展开求得F-1的近似表达式，往往因计算工作量过大而无法实现。

函数变换法一、原理若随机变量X具有密度函数f(x)，分布函数为F(x)，随机变量Y是X的函数Y＝g(X)，当逆函数x=g-1(y)=h(y)存在且连续的一阶导数时，随机变量Y的密度函数为p(y)=f[h(y)]|h’(y)|，由此算出其分布函数G(y)＝F[g-1(y)]，再用逆变法可得抽样公式Y=g(X)。函数变换法是逆变法的推广。函数变换法生成随机变量的一般步骤是：①产生独立的F(x)的随机数x1，x2，…xn；②由抽样公式，得随机数列yi=g(xi)（i=1，2，…，n）。二、正态分布正态分布的密度函数为分布函数为通过Z变换，有令则有无法用逆变法直接从上式求得随机变量Z。可采用Box-Muller的函数变换法。先产生两个独立的标准正态分布Z1和Z2，因为Z1、Z2独立，故它们的联合密度函数为令代入，则有联合密度函数f(r,θ)为r的函数与θ的函数之积,即

f(r,θ)=f1(θ)•f2(r)其中：只要产生随机数u1，u2∈U(0,1)，使θ=2πu1，即得通过Z变换

得到一般正态分布N(μ，σ2)随机变量：

X1＝μ＋σZ1

X2＝μ＋σZ2三、对数正态分布若X～N(μ，σ2)，则称Y=eX所服从的分布为对数正态分布，记为LN(μ，σ2)。其均值和方差是：利用变换公式Y=eX，则可由N(μ，σ2)的随机数生成LN(μ，σ2)随机数。具体步骤是：

①产生N(0，1)分布随机数ui；

②计算xi

=μ+σui；

③令

，则得对数正态分布LN(μ，σ2)随机变量。用函数变换法由N(3，0.12)分布的随机数生成对数正态分布变量。解：用Excel的随机数发生器产生N(0，1)分布随机数ui，然后由xi

=3

+0.1ui得N(3，0.12)分布随机数，再据

产生对数正态分布随机变量，结果如下。生成随机变量的方法主要有：逆变法、函数变换法、舍取法、卷积法、组合法等。常用分布的随机变量可在仿真软件包中直接找到，在实际使用时调用相应的函数即可。仿真输入数据分析输入数据的正确性直接影响仿真的结果，正确收集和分析数据是系统仿真的重要组成部分。输入数据分析的一般步骤是：1、数据收集；2、分布识别；3、参数估计；4、假设检验。数据收集足够数量的有效原始数据，是仿真有效的基础。在收集过程中，当数据不足时，应补充收集；而与仿真无关的数据则应予剔除。离散随机系统仿真以统计方法为基础，而统计方法要求样本独立地取自某个总体，因此在分析仿真输入数据时，首先要确定这些样本是否独立。判断样本独立性的方法可有：1、参数方法；2、非参数方法。独立性检验一、散布图法利用样本点的散布图，可以粗略地判断样本是否独立。如果随机变量相互独立，其样本散点图随机地分布在座标平面1030507090-300-200-1000100200yy

散布图二、相关图法相关图法是依据样本相关系数

的图形判断样本是否独立。若样本X1，X2，……，Xn独立，则相关系数

若样本相关系数

不等于零，则样本X1，X2，……，Xn不独立分布的识别随机变量的分布类型，可利用先验知识，关于所研究系统的实际了解，输入数据的物理意义，假设随机变量分布的类别。但在实际工作中，更常遇到的数据集合是无法用先验信息确定输入的随机变量的分布，必须利用其它方法。本节介绍：1、直方图方法；2、Q-Q图法；3、数字特征方法。一、直方图法直方图法是根据所收集的有效数据，绘制直方图，然后分别与几类不同的概率密度函数曲线比较，确定直方图与哪一类分布的密度函数接近或相似。例如：某银行为提高服务窗口的工作效率，决定用仿真研究窗口的设置。为此，需确定在固定时间内顾客到达间隔时间X(min)的分布形式。在固定的时间内有221名顾客来到，每相邻的两个顾客到达的时间间隔，从小到大排列如表顾客到达的时间间隔二、分位点-分位点（q-q）图q-q图建立在yi是X的（j-1/2）/n分位点的估计，即

的近似的基础上。若选定的分布F是X的一个合适的分布，那么yi关于

的图形接近一条直线，如果F的参数值合适，该直线的斜率为1。汽车装配线上车门的安装时间如下序号安装时间序号安装时间序号安装时间序号安装时间199.79699.5611100.1716100.332100.267100.411299.981799.833100.238100.2713100.0218100.47499.55999.621499.651999.82599.961099.915100.062099.85

q-q图用Excel作q-q图虽然可省去繁杂的计算，但还是比较麻烦。现在已有许多数据软件，例如@risk、SPSS等软件包可方便地用来产生q-q曲线。用@risk画q-q图三、数字特征法

常用数字特征部分地描述分布的一些特征，而这些数字特征又可用样本数据来估计，因此可用数字特征概略地估计输入数据的概率分布。常用数字特征及其点估计一些常用分布的莱克西斯比率和偏差系数参数估计为在仿真过程中进行抽样，确认输入数据分布类别后，还必须利用观测的样本数据估计参数的值。参数估计有许多不同的方法，以下介绍常用的极大似然估计法。极大似然法的步骤如下：1、写出似然函数

2、求出lnL及似然方程

3、解似然方程求出参数θj的估计值

。单位时间打入电话总机的电话次数是服从泊松分布的随机变量其中xi是单位时间打入的电话次数，λ是未知参数，其似然函数为或于是有则因此，λ的极大似然估计为一分钟内打入电话数均值=4.895假设检验对输入数据进行独立性判别、分布形式的识别、参数估计以及分布函数或密度函数的估计之后，还需要用所得样本的信息对分布形式、所估计的参数等进行假设检验，判断假设的正确性。1、柯尔莫哥洛夫—斯米尔诺夫（Kolmogorov-Smirnov）检验；2、分布拟合的χ2

检验一、柯尔莫哥洛夫—斯米尔诺夫（K-S）检验K-S检验是在每一点上考察样本经验分布函数与总体分布函数间的偏差，适用于理论分布F0完全已知的连续分布函数情况。具体步骤是：1、假设2、随机数x1，x2，……，xn按大小排列，记为x（1），x（2），……，x（n）；3、计算K-S检验的统计量

的值

其中4、对给定的水平α，查K-S分布表的临界值Dα(n)：

5、若

则接受假设H0，即认为观测数据与理论分布拟合较好；否则，拒绝假设，认为观测数据与理论分布不相符。

用柯尔莫哥洛夫—斯米尔诺夫检验法检验如下一组容量为10的样本：是否接受分布函数为：

均匀分布的假设X1X2X3X4X5X6X7X8X9X100.6210.5030.2030.4770.7100.5810.3290.4800.5540.382计算K-S统计量临界值为0.409，Dn＝0.29<0.409，接受零假设，认为服从分布函数如上的均匀分布。服从均匀分布二、分布拟合的χ2

检验K-S检验考察每点处的经验分布与理论分布的差异，χ2

检验分组处理样本区间，用K.Pearson提出的统计量χ2来检验总体是否服从一个预先给定的连续或离散的分布。χ2

检验的具体步骤是：1、假设

2、将实数轴分成若干互不相交的区间[a0,a1)，[a1,a2)，……，[ak-1,ak)

3、计算概率pi=P{ai-1<X≤ai}=F0(ai)－F0(ai－1)，（i=1,2,……,k），由此计算理论频数npi；

4、计算落在每个区间[ai-1,ai)的样本观测数，即经验频数fi；

5、计算皮尔逊统计量

6、由χ2分布表查给定水平α的临界值

；

7、比较χ2与临界值

，若

则接受零假定，反之则拒绝零假设。用χ2检验银行服务窗口的数据是否符合指数分布？因为λ＝2.380695，理论分布函数为将区间[0,∞）分为如下10个子区间

区间0-0.0410.041-0.080.08-0.120.12-0.200.20-0.27经验频数2422222421概率0.09300.08040.07510.13030.0954理论频数20.459117.692716.518028.671120.9776统计量χ20.61281.04861.81930.76100.0000区间0.27-0.390.39-0.520.52-0.670.67-0.9750.975-∞经验频数2421192320概率0.13070.10520.08710.10470.0982理论频数28.746823.140219.157523.042221.59471、计算各区间的概率和理论频数

2、计算到达间隔的统计量χ2

3、由表查临界值

χ2＝5.3427

<＝15.507接服从λ＝2.380的负指数分布的假设。离散事件系统仿真采用的MonteCarlo方法，其本质是统计抽样试验，仿真过程中输入数据是来自总体的样本，样本的有效性决定决定仿真结果的有效性，不可靠的输入数据会导致错误的建议。因此，必须对输入的数据进行分析，输入数据分析的内容是：数据收集，分布识别，参数估计，假设检验。

Spreadsheet建模

第一节Excel基础Excel是MicrosoftOffice的一个组件，界面友好、功能强大，应用广泛。本节介绍与建模有关的Excel2007的知识。公式与函数spreadsheet建模涉及各种公式，公式是对工作表中的值计算的等式，由函数、单元格、常量以及运算符四个部分组成，如=sin(D2)^2一个公式可包含以上四部分中的部分或全部。

Excel的内置函数Excel内置有关财务、数学与三角函数、统计、查找与引用、数据库、文本、逻辑、信息、工程等类型约二百多个函数。每个函数均有一个唯一的名称，引用函数格式是：

=函数名（）输入的方式有三种：1、利用插入函数对话框；2、利用屏幕提示3、直接输入

在创建公式时，常会用到嵌套函数，在Excel2007中，最多可以嵌套64层函数。例.在下表中有随机变量x1，x2，…，x7

，每个变量样本量为5，自x2起xi6的取值按以下规则：当前一组变量的5个样本的均值大于30，xi6为本组样本均值，若前一组变量的5个样本的均值小于30而大于25，则为前一组均值，否则为0。x1x2x3x4x5x6x721.8439335.2037133.9512330.1651127.3445838.7633924.3757431.7123935.8421639.0484330.6015238.5485438.5174122.2919424.6028032.7317125.2607820.6274634.5597725.979822.1381322.1973327.0290235.662135.9221220.0543232.7500229.8849531.3211522.9700626.7012529.8147536.5306621.2921523.0488要创建公式“只有当前一组数字(F2:F5)的平均值大于30时，才计算本组数字(G2:G5)的均值,若这组数字(F2:F5)的平均值小于30而大于25，则为前一组之均值，否则为0”，在编辑栏输入公式：

=IF(AVERAGE(A2:A5)>30,AVERAGE(B2:B5),IF(AVERAGE(A2:A5)>25，AVERAGE(A2:A5)，0))数组与模拟运算（一）数组（一）数组数组是以一个范围内的整个数组作为运算对象的特殊公式。（二）模拟运算表将单元格看成是其它单元格的函数，在单元格实现复杂的规划求解，甚至是模型的运算这样一种一功能的工具。

数组的形式比较特殊，是在公式的等号“=”之前与公式结尾处加一组花括号“{}”，这个花括号是在选取的范围的条件下，正确输入公式后，按组合键Ctrl+Shift+Enter后自动生成，而不能在键盘上输入的。数组中的任一部分不能单独修改。利用Del键可以删去数组。若不慎删去数组的一部分，可用组合键Ctrl+Shift+Enter删去数组。

数组例：将下表中前6行的数据精确到小数点2位。45.9788849.7802733.8816533.5804347.910741.3296936.792233.4046449.0862846.9374749.3499637.2637133.9063740.2850424.8341337.3268834.8493932.0093442.1573549.4140441.7929645.456137.3296322.3401622.7347628.281227.4123439.823626.8382827.0854825.3825539.6258425.7231437.245429.3551429.6453733.9805335.16834.5976125.6123544.7996525.1216231.7804530.1068121.689237.8469229.8843327.5624922.19825

在Excel2007中，模拟运算表称为数据表。数据表作为一个函数，虽然可有多个自变量，由于表格是二维的，即不同的输入只能通过行、列来实现，因此可控制或调节的自变量至多有二个，即只能有一维和二维数据表。模拟运算表月

份0123456股票A价格25.9928.5519.0321.9215.9741.2634.67股票B价格47.0672.8863.6123.2147.7160.8263.26月

份789101112股票A价格27.7929.5729.629.8729.4935.17股票B价格60.547.243.4864.8348.8653.24投资组合的标准差和期望受益率曲线。有A、B两只股票过去12个月的价格如下:计算并画出不同组合的标准差收益率曲线。一维数据表的例连续收益率公式

月期望收益为12个月的收益率平均值，年期望收益等于12乘月期望收益，按公式

计算组合收益率在H3：J24选定区域二维数据表的创建与一维数据表类似，不同的仅仅是要输入“输入引用列的单元格”和“输入引用行的单元格”。二维数据表利用数据表编制变量精确到小数点后二位的正态分布表。审核

审核是电子表格建模的基本要求之一。Excel的公式审核功能提供了有效、便利的检查公式的工具——“切换（toggle）”。公式审核还提供用图形直观地追溯单元格之间关系的方法，称为追踪单元格。追踪单元格可分为：追踪从属单元格和追踪引用单元格两种。

标注与域名在电子表格模型中，单元格或区域使用域名或标注，可使公式和模型容易读懂和理解。

定义域名的方法如下，在选定单元格或区域后，单击【公式】、【定义名称】，然后在对话框【新建名称】中输入域名，【确定】.宏在实际问题中，常会遇到通过改变模型的一些参数，反复运行，比较不同的方案，寻求最满意的结果。Excel提供只要单击一个按钮即可自动按建模者意愿运行的方法，这就是用VAB创建“宏”。连锁快餐店有12家分店，分发年终奖金的依据是用投入、产出表示的运营效率。有关12家分店的运营数据如下：分店123456789101112利润5.987.184.975.323.394.952.896.46.016.945.868.35满意度7.79.79.37.77.87.98.69.17.38.88.29.6保洁度9299988794889010089899397工作小时4.746.385.043.613.455.252.367.096.497.365.466.58运营成本6.757.426.356.344.436.313.238.697.289.076.698.75快餐连锁店效益数据

分店的效率指数为分店的产出加权平均值与投入加权平均值之比，即：

其中Oij为第i分店第j种产出的数值，Iij为第i分店第j种投入的数值，uj、vj分别是第j种产出和第j种投入的权重，hi为第i分店的效率指数。≤1

投入的加权平均值等于1，即作如下变换

,,问题转换如下线性规划问题

对第1分店的效率，应解线性规划

第二节Spreadsheet建模Spreadsheet——“电子制表软件”或“电子数据表”。

Spreadsheet建模的优点1、简单易用2、功能强大3、易于理解建立电子表格模型的一般步骤设定目标；手工试算；建立表格

建立模型

测试、分析模型逻辑关系评估解和/或用solver优化模型需要扩充、完整？模型是否完善？否或否是建模完成是1、设计。就是设定电子表格的结构。可进一步分为确定目标、试算以及草拟电子表格三部分。2、建模。考虑好电子表格的合理编排后，就可以打开Excel工作表，建立模型。3、测试。通过在可变单元格中输入数值，看输出结果是否与所知道的输出单元格的数值一致。4、分析。运行Excel的Solver(在工具的下拉式菜单中的规划求解)求解及优化。年份12345678910净现金流-6-2-4243-46-18例：某公司计划部门对今后10年的现金流预测如表

公司净现金流预测数据单位：百万数据显示公司在近三年收入不足维持营运及成本，几年后可转入盈利。目前，公司有100万现金余额，为确保公司正常营运，至少需要保持50万的现金储备。公司有业务往来的银行可提供利率为5%的10年后还本的长期贷款，以及利率为9%的一年期的短期贷款。公司应采取怎样的策略，既满足运行需要，又使支付利息最少？

解：（一）设计公司的目标是筹集资金，支付最少的利息。

诸变量之间最重要的关系，是期末余额与期初余额、现金流、贷款额、利息等之间的平衡：

期末余额=期初余额+现金流+贷款-利息支出-贷款偿还第i年的期末余额与期初余额、现金流、贷款额、利息等之间的平衡关系可写成：

根据以上分析，草拟表格如下：年份现金流贷款支付利息还本期末余额最小现金余额长期短期长期短期长期短期1-64210.52-254

0.052

0.0921.620.53-44

0.055

0.0950.5424

0.050.5544

0.050.5634

0.050.57-44

0.050.5864

0.050.59-14

0.050.51084

0.050.5114

0.0540.5建模先建前两年的小型模型建立电子表格模型应注意的事项1．数据结构精心编排2．公式数据应分开输入3．输入的公式应简单明了4．不同功能的区域应标识清楚电子表格模型的调试

1．数值试算。2．检测域名。3．审核公式。

4．追踪单元格。Spreadsheet建模的例一、目标规划二、财务模型三、预测模型四、判别分析

超市配送中心S1S2S3S4供应量D1583410D274524D3269116需求量1210612华盛集团在某地有四个大型超市，由三个配送中心配送货物。已知各配送中心与超市的供需数量及运价如下：表3.3.1供需数量及运价表公司的经营目标是：

1、所有配送中心的货物应全部运出：

2、每个超市至少得到需求量的80%；3、超市S3的需求必须满足；4、由于道路修建，应尽量减少D1到S2的货运量；5、使总运费最小。请研究设计满足以上目标的方案。Min()

公司财务计划某公司准备用仿真制定今年的财务计划。根据财务部门分析销售成本（CGS）等于销售额（Sales)的85%，预计销售成本不少于上年销售成本的90%。需支付的利息为负债额的9%。为保证公司营运，应保留盈余（上年的保留盈余与纯利之和），现金（上年的现金与纯利之和再加债款增量）不低于100000元。公司固定资产为500000元。已知上年的销售额为1100000元，债款为200000元，预计今年销售额为1400000元，债款为200000元。作为一个例子，先研究一个简化的企业财务模型。假定企业的财务关系是：（1）销售成本等于销售额的8%，但不低于上一年销售成本的90%；（2）毛利等于销售额减销售成本；（3）利息等于债款的9%；（4）税前利等于毛利减利息；（5）所得税等与税前利的48%；（6）纯利等于税前利减所得税；（7）收益留成等于上一年度的收益留成加本年度纯利；（8）现金等于上一年度现金加本年度纯利和本年债款增量，；（9）总资产等于现金加固定资产，假定固定资产每年均为500000；（10）总负债等于债款加收益留成。

据有关部门统计预测，今年销售额为1000000，今后三年的销售额分别为，700000，1100000，1400000。根据以上关系，建立电子表格模型财务模型成本率85%现金

100,000利率9%固定资产

500,000税率48%年份2010201120122013损益计算销售额

1,000,000

700,000

1,100,000

1,400,000销售成本

850,000

765,000

935,000

1,190,000毛利

150,000

-65,000

165,000

210,000利息

18,000

22,075

15,058

5,487税前利润

132,000

-87,075

149,942

204,513税额

63,360

-41,796

71,972

98,166纯利

68,640

-45,279

77,970

106,347资产平衡现金

100,000

100,000

100,000

100,000固定资产

500,000

500,000

500,000

500,000总资产

600,000

600,000

600,000

600,000负债

200,000

245,279

167,309

60,962留成收益

400,000

354,721

432,691

539,038总负债

600,000

600,000

600,000

600,000三、三次指数平滑预测模型三次指数平滑的目的与二次平滑一样，是建立预测模型平滑。它是一种非线性平滑，可更好地预测拐点。例：通达公司为做好财务决策需要预测销售额，已知某产品2004年以来的季度销售额（表3.3.2），试用三次指数平滑为通达公司预测下一年度销售额。表

通达公司产品销售额单位：百万年度200420052006季度123412341234销售额240.3222.8243.1222.2220.6218.7234.5248.6261.0275.3269.4291.2年度200720082009季度123412341234销售额301.9258.5286.6260.5298.5291.8267.3227.9303.5313.3327.6338.3布朗三次指数平滑的预测模型为：

其中α是平滑指数，

分别是一次、二次和三次平滑值，模型中的参数at，bt，ct由以下公式来确定

四、判别分析从本质说,判别分析是对离散变量进行分组预测的统计技术。可进一步分为对于两个组分类和用多组分类的判别分析法（MultipleDiscriminantAnalysis）。例:ACME公司每年要对员工进行“机械智商”和“口头表达能力”的测试，有关测评数据如表:序号12345678910类别1111111111机械智商44.743.242.540.241.939.338.638.137.936.8口头表达36.142.030.842.637.133.840.328.544.742.6序号11121314151617181920类别1222222222机械智商35.839.037.536.135.634.432.831.429.530.1口头表达35.233.136.430.827.239.035.229.437.933.9表中类别1为合格，2为不合格。假设机械智商及口头表达能力与合格性之间的关系线性关系，即参数b0、b1、b2可用最小二乘法估计。机械智商及口头表达能力与合格性之间的关系为：若公司新来求职者的测评数据如表新员工编号123456机械智商40.541.636.443.829.536.3口头表达33.440.432.928.734.427.4新求职者测评数据表蒙特卡罗方法与随机数

蒙特卡罗（MonteCarlo）方法，又称随机抽样法或统计实验法，是基于对概率模型的观察或抽样试验得到问题近似解的方法。本章讨论：

1、蒙特卡罗方法的基本思想；

2、随机数的产生和检验。

蒙特卡罗（Monte-Carlo）方法

蒙特卡罗方法以概率统计理论为其主要理论基础，以随机抽样为主要手段。先建立一个概率模型或随机过程，使它的参数等于问题的解；然后通过对模型或过程的观察或抽样试验来计算所求随机参数计算所求随机参数的统计特征，最后给出所求解的近似值。

若射击弹着点到靶心的距离r，则得分为Y，Y与r的关系为g(r)，即

r用统计试验的方法求其数学期望

N次射击的弹着点股票价格的预测研究

（ITO）过程

维纳过程（Wienerprocess）

服从正态分布

σ＝0.03，Δt=0.01,股票的初始值S=30元，用变换

每周期开始价随机样本v1∈N(0,1)随机样本v2股价变化ΔS每周期结束价30.000-0.12675-0.0026-0.0780829.92229.9220.7216090.0228480.68366430.60630.606-1.08988-0.0315-0.9639729.64229.642-0.60426-0.01693-0.5017729.14029.1401.9316120.0591481.72357530.86330.863-0.12521-0.00256-0.078930.78530.785-1.266-0.03678-1.1322629.65229.6521.5679780.0482391.43040631.08331.083-0.26258-0.00668-0.2075630.87530.8750.0898150.0038940.12024230.99530.9951.0504160.0327121.01393632.009（1）构造或描述概率过程。（2）实现抽样。（3）建立估计量。实施蒙特卡罗仿真的主要步骤：蒙特卡罗仿真的优点：（1）收敛速度与问题的维数无关。（2）受问题的条件限制的影响小。（3）程序结构简单、清晰，便于编制和调试。（4）对于一些物理问题，具有其它数值方法不能替代的作用。随机数的生成和检验

随机数就是按随机方法而生成的数码。即0，1，2……9这十个数字出现的机会等概率，排列的顺序随机。均匀分布随机数随机数的产生方法0123468932765140675454324678568087586545均匀分布随机数f(x)F(x)xx1010111

均匀分布概率密度及分布函数

随机抽取等概且服从均匀分布的小数（从0到1，但不包括1）称为均匀随机数

产生均匀分布随机数方法

手工方法2.随机数表

物理方法4.数学方法

手工方法掷骰子、抽签、发纸牌或从“搅拌均匀的容器”中模编号球的方法等，但效率很低。随机数表物理方法对随机数表法的一种改进办法是在计算机上装一台物理随机数发生器，把具有随机性质的物理过程（如以放射性物质为随机源），在计算机上直接转换成随机数字。这种方法虽然产生了随机数，但无可追溯性。由此得到的随机数一去不返，不能重复出现，因此无法再用原来的随机数试算或检查。又因为对设备的要求较高，而实用价值较低

数学方法利用数学递推公式，在计算机上产生随机数，是目前使用较广、发展较快的一种方法。由于这种方法产生的数，只能近似地具备随机性质，因此称为伪随机数（pseudorandomnumber）。产生伪随机数的方法

（一）历史上曾经采用过的方法1.平方取中法(mid-squaremethod)

平方取中法计算简单，但有明显得缺点，一是无法说明用什么样的种子数可保证有足够长的周期，二是容易退化到某一常数或零，而一旦有一个零，则以后的数全为零以21为初始种子数，用平方取中法求随机数。xixi2xi+1ui21449364986300441193686494096008100640036000900004493649863000.440.930.640.090.080.060.030.000.002.乘法取中法(mid-productmethod）

设初始种子为12和34，用乘法取中法求随机数。xixi•xi+1xi+1ui1234403644585519472100408136014491584255231901045007600280014000200004036445855194721000.400.360.440.580.550.190.040.070.020.010.000.003.常数乘子法（constantmultipliermethod）

常数乘子法的周期较长，均匀分布特性较好，但仍有平方取中和乘法取中的缺点——容易退化。同时，成功与否取决于所选常数。xiki•xixi+1ui341870870.87874785780.78784290290.29291595590.59593245240.24241320320.32321760760.76764180180.1818990990.99995445440.44442420420.42422310310.31311705700.07703850850.85854675670.67673685680.68683740740.7474407070.07341870870.87874785780.78选择常数k=55，初始种子数为34，用常数乘子法求随机数。

（二）线性同余法1.混合线性同余法

2.乘法线性同余法3.二次同余法

4.加法同余法

用线性同余、混合线性同余、二次同余及加法同余，计算随机数。其中：a=21，b=59，c=11；模m=128；种子数是：x0

=34,x-1

=51。参数m，a，c满足条件：①m与c互素；②如果q是m的一个素因子，则q也是a-l的因子；③如果m能被4整除，则a-1也能被4整除。通常取m=2b,c为奇数，而a-1可被4整除，可得到满周期。随机数的检验

由于计算机产生的随机数，是按某一迭代公式计算出来的，因此实际上由此而得的数列是确定性而不是随机的。为此，通常要对随机数进行测试，观察它与均匀分布U(0,1)的类似程度。随机数检验可分为:参数检验均匀性检验独立性检验

检验随机数的均值、二阶矩及方差等估计量与均匀分布U(0,1)的理论值有无显著差别

参数检验用参数检验由某随机数发生器生成的100个随机数的均匀性。

0.38200.10070.59650.89910.88460.95850.01450.40740.86320.13860.24500.04550.03240.16410.21960.01710.28500.34310.55360.35740.37180.35560.91030.46600.42620.30390.97570.80670.99120.25630.95170.05340.70500.81650.97250.46630.30020.75020.35150.77570.07430.19840.06410.35830.48700.51120.37350.98590.04070.23070.00500.92610.10030.25670.77570.67960.80910.72430.08510.13230.75620.62650.17370.40480.55230.71150.55520.18120.97030.68690.52880.79670.80570.26220.17800.86680.11480.05950.76160.73840.98630.92560.90390.54500.50080.67500.48980.14580.03800.79630.67160.73170.58450.15220.89220.37780.20050.20580.33400.3251

检验随机数在各子区间的频率与理论值N/n有无显著差异。常用的检验方法有:频率检验（拟合优度检验）

累积频率检验（柯尔莫哥洛夫－斯米尔诺夫检验）

均匀性检验频率检验是按以下方法，对随机数独立且服从U(0,1)分布的假设作显著性检验。具体步骤是：①将区间[0,1]分成k个不相交的等长子区间②计算经验频数③计算统计量④判断统计量χ2是否小于临界值在5％的显著性水平上，用χ

2检验所给的100个随机数的均匀性。将区间[0,1]分成10个等间距的子区间，各子区间的随机数个数ni

如表4.2.5:

经验频数表子区间i12345678910

ni121291379512912据表的经验频数计算统计量χ2

χ2＝6.2<16.92

累积频率检验又称柯尔莫哥洛夫－斯米尔诺夫检验，是连续分布的拟合性检验（简称K-S检验）。按以下步骤对样本的经验分布函数Fn(x)与总体分布函数F(x)进行显著性检验。

①将随机数按大小顺序排列成u(1)，u(2)，…，u(n)，并计算随机数的经验分布函数Fn(x)②计算

③计算统计量④根据显著性水平α及样本大小n查K-S表检索临界值⑤比较统计量D及临界值

【例4.2.7】在5％的显著性水平，用K-S法检验例4.2.5前10个随机数的均匀性。解：将随机数按大小顺序排列，得相关数据如表

Dn+=0.161,D-n=0.263，

Dn=max{0.161,0.263}=0.263

D0.05(10)=0.41

i12345678910u(i)0.0150.1010.1390.3820.4070.5970.8630.8850.8990.959i/n0.10.20.30.40.50.60.70.80.91.0i/n-u(i)0.0850.0990.1610.0180.0930.003-0.163-0.0850.0010.041u(i)-(i-1)/n0.0150.001-0.0610.0820.0070.0970.2630.1850.0990.059

独立性检验就是检查一个序列的随机数之间是否存在相关性。常用检验的方法有：自相关检验

组合规律检验

连贯性检验

独立性检验

利用自相关系数值的大小来判断随机数的独立性，因此也称相关系数检验。具体步骤是：

①计算j阶自相关系数的估计值

②构造统计量vj（j=1，2，…，m）

③根据显著性水平α，比较zα/2和统计量判断待检随机数是否具有独立性。

>zα/则拒绝相关系数ρj=0的假设

自相关检验

通常用扑克检验（pokertest）检验样本中重复出现数字的频率，从而检查其组合规律性。

组合规律检验

若Ak为随机数中有且仅有k个数字相同的事件，Ak出现的频数为mk。则三种可能事件发生的理论频数为μk=nPk。构造统计量

由于统计量χ2渐近服从χ2(2)分布。因此可利用χ2分布检验随机数列的组合规律性来判断其独立性。å=-=3122)(kkkkmmmc有1000个用三位小数表示的随机数，经过统计知有680个三位数都不相同，289个恰有两位数字相同，31个三位数字都相同。试检验这些数据是否相互独立。解：三种不同组合情况的有关数据如下：据此计算统计量,若显著性水平α＝0.05，查分布表得，因为＝47.66>5.991，拒绝符合组合规律，表明这1000个随机数不独立。类型概率观测频数mk

理论频数μk

A10.7268072016002.22A20.272892703611.34A30.01311044144.1连贯性检验通过考察数字的排列来检验随机数独立性假设，即检验随机数列先后出现的随机性，重点是检验其连贯现象是否异常，从而判断这些数据是否独立。连贯检验可分为：1、升降连检验2、正负连检验。升降连检验由u1，u2，…，un，构造数列{vi}，vi＝ui-ui-1(i=2,3,…,n)，根据数列{vi}的符号，分成升降两类连，vi连续为正的连称为升连，vi连续为的连称为降连。用连来表示随机数增减及其长度的变化规律的方法称为升降连检验。ui0.7440.9060.2640.2310.2570.0070.1590.3630.7140.070Vi0.162-0.642-0.0330.026-0.2500.1520.2040.351-0.644+--+-+++-ui0.3690.2680.8980.13620.0860.4730.7080.2450.4380.132vi0.299-0.1010.630-0.762-0.0500.3870.236-0.4630.193-0.306+-+--++-+-ui0.4420.2270.49390.5810.6240.5980.7450.8060.6180.703vi0.310-0.2150.2670.0870.043-0.0270.1480.061-0.1880.085+-+++-++-+ui0.2740.5580.47430.1630.4350.9830.1970.6360.6400.530vi-0.4290.284-0.084-0.3110.2720.548-0.7860.4390.004-0.110-+--++-++-ui0.4390.4090.94630.9740.7160.0640.7910.7110.5520.496vi-0.091-0.0300.5380.027-0.258-0.6520.727-0.080-0.159-0.056--++--+---{vi}及升降连升连、降连共32个,其中长度为1的有19个,长度为2的有9个,长度为3的有4个T=32，n=50，统计量

又zα/2＝z0.025=1.96，可见|Z|<z0.025，所以不能拒绝独立性假设。正负连检验正负连是将数值在均值以上的标以“＋”，表示正连；在均值以下的标以“－”，表示负连。然后，用类似升降连检验的方法进行检验。

假设n1和n2分别是在均值以上和以下的数据个数，T为连的总数，其均值和方差分别是

，

（5.3.11）当n1或n2大于20时，T近似服从正态分布，统计量

（5.3.12）服从标准正态分布N(0,1)。同样地，给定显著性水平α，若|Z|≤zα/2，则不拒绝独立性的假设。

显著性水平α＝5％，用正负连法检验上例所给的随机数列的独立性。ui0.7440.9060.2640.2310.2570.0070.1580.3630.7140.070vi0.2580.420-0.222-0.255-0.229-0.479-0.327-0.1230.228-0.416符号++------++ui0.3690.2680.8980.1360.0860.4730.7080.2450.4380.132vi-0.117-0.2180.412-0.350-0.400-0.0130.222-0.241-0.048-0.354符号--+---+---ui0.4420.2270.4940.5810.6240.5970.7450.8060.6170.703vi-0.043-0.2590.0080.0950.1380.1120.2590.3200.1320.217符号--++++++++ui0.2740.5580.4740.1630.4350.9830.1970.6360.6390.529vi-0.2120.072-0.011-0.323-0.0510.497-0.2890.1500.1540.044符号-+---+-+++ui0.4390.4080.9460.9740.7160.0640.7910.7110.5520.496vi-0.047-0.0770.4610.4880.230-0.4220.3050.2250.0660.011符号--+++-++++n1=25，n2=25，T=19，并有

因为|Z|=1.57175<1.96，故不拒绝独立性假设。

Monte-Carlo方法是以概率统计理论为理论基础，随机抽样为主要手段，程序结构简单、清晰，便于编制和调试，对于一些物理问题，它具有其它数值方法不能替代的作用。因而在许多领域得到广泛的应用。随机抽样关键在于生成随机数，随机数就是[0，1]区间均匀分布的随机变量。随机数可用手工方法、随机数表、物理方法和数学方法来产生。

用数学方法来产生随机数，通常借助计算机由公式生成，因此是伪随机数。在仿真时只需直接调用仿真软件的随机数函数即可得到随机数。用公式产生的伪随机数，必须要经过测试，检验它的均匀性检验及独立性检验，是否符合参数分布。运用

Spreadsheet仿真Spreadsheet仿真模型是电子表格模型的一种。它与其它电子表格模型的区别在于：它明确地将不确定性结合到一个或多个输入变量。在模型运行时，可对这些变量随机输入不同的值，观测模型如何岁输入的变化而变化，回答各种假设性分析问题。Spreadsheet仿真模型

如果例3.1.1未来10年的现金流的预测值不是确定值，而是服从正态分布的随机变量，其均值为表8.1.1提供的数据，标准差为均值的5%。

公司净现金流分布

单位：百万年份12345678910均值-6-2-4243-46-18标准差0.30.10.20.10.20.150.20.30.050.4数据驱动仿真模型Spreadsheet仿真不但需要建立电子表格模型，还必需使输入的随机变量的概率分布符合实际。因此必须解决以下二个关键问题：一、在所讨论的问题选择某种分布的理由是什么？二、仿真模型中可使用的分布有哪些，如何在模型中调用？PortaCom公司设计了一种外观新颖、打印质量高的新型便携式打印机，预计能在市场占有很高的份额。经对市场和财务的初步分析，零售价、第一年的行政管理费和广告费为：

零售价＝249（美元）

行政管理费用＝400000（美元）

广告费用＝600000（美元）

由于劳动力费用（工资）、零件材料费用及第一年的需求量无法确切得知，根据预测每台工资为43～47美元，每台零件材料费为80～100美元，第一年的需求在5000～15000台之间。因为公司银根较紧，对亏损的可能性极为关注，因此要求进行利润分析以便决策。ProtaCom仿真流程单元格输入公式作

用C17=NORMINV(RAND(),$G$10,$H$10)预测需求D17=VLOOKUP(rand(),$B$10:$D$14,3)预测工资E17=$G$13+($H$13-$G$13)*rand()计算材料费F17=$C$5行政管理费用G17=$C$6广告费H17=B17*C17销售总额I17=D17+E17+F17+H17总成本J17=H17-I17利润@RISK定义分布（DefineDistributions)添加或编辑概率分布函数输出（AddOutput)添加仿真输出插入函数（InsertFunction）插入@RISK函数相关性（DefineCorrelations）确定随机变量相关性分布拟合（DistributionFitting)拟合数据分布分布画家（DistributionArtist）手绘分布曲线模型视窗（ModelWindow）显示输出单元格及分布函数——迭代次数复选框，给定每次仿真运算次数；——仿真次数复选框，确定仿真的次数。

仿真设置迭代次数、仿真次数、抽样类型等仿真设置重算功能等同于Excel的重算键F9图形输出自动显示输出图形仿真概要结果自动显示概要结果显示模型运行时展示或隐去模型在线更新运行时即时更新@RISK窗口运行仿真单击图标则开始运行仿真模型结果高级分析点击右下角按钮

，可从下拉菜单可选择目标搜寻、应力分析或高级敏感性分析。浏览结果,显示与所选单元格有关的结果图形输出结果筛选结果统计汇总数据敏感性分析情景分析结果的Excel报告包括数据表、图形等切换@RISK库打开设置显示平行窗口，仿真文件，清除@RISK数据，卸载@RISK@RISK帮助

输入随机变量@RISK有三种输入随机变量的方法：1、直接利用DefineDistributions2、利用InsertFunction3、直接输入@RISK函数用DefineDistributions输入利用InsertFunction直接输入@RISK函数@RISK函数均冠有前缀Risk，在@Risk5.5中，一般分布函数（common）16个，进而分为离散8个、连续33个、特殊5个、参数备用34个（递增和递减各17个）。输入@RISK函数的方法有：1、单元格输入公式2、单击【公式】选项卡、插入函数图标

，用【插入函数】对话框选择函数3、利用Excel2007提供的提示功能。@RISK函数1、用RiskNormal产生正态分布随机变量2、用RiskUniform产生均匀分布随机变量3、用RiskTriamg产生三角分布随机变量4、用RiskBionom产生二项分布随机变量5、用RiskBionom产生二项分布随机变量6、用RiskSimTable函数提供多个随机变量常用分布的@RISK分布函数

仿真设置，是指确定迭代次数（重复仿真的次数）、仿真次数以及取样方式。仿真设置General表View表Sampling表Macros与Convergence表

仿真运行结束时，@Risk5提供统计报表、仿真数据、敏感性分析以及情景分析帮助解释仿真结果。工具棒

表示以上四种报告的图标，点击其中的一个图标，即可显示相应报告的窗口。报告生成

光华书店须在9月份确定订购挂历数量，按历年统计，挂历的需求量为服从均值为200，标准差为50千本的正态分布的随机变量，供应商能提供的挂历是125、200及250千本的三角分布的随机变量。若供应商不能满足书店订购数，每本的价格为7.5元，不然每本的价格为7元，挂历的销售价格为每本15元。在1月15日后，挂历的销售价格为每本5元，若挂历到3月1日还卖不出，只能作废纸处理。对剩余挂历的潜在需求量为剩余挂历的0%、75%间均匀分布的随机变量。书店可采取的订购策略有150、175、200、225、250千本五种。用仿真的方法评价各种订购策略。@RISK建模单元格

输入公式公式含义A30=Roun