Monte-Carlo模拟

1、蒙特卡罗方法蒙特卡罗方法（Monte Carlo Simulation）曼哈顿计划曼哈顿计划Buffon投针实验投针实验大数定律大数定律基本思想：当所求问题是某种随机事件出现的基本思想：当所求问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过概率，或者是某个随机变量的期望值时，通过某种某种“试验试验”的方法，以这种事件出现的频率的方法，以这种事件出现的频率估计估计该该随机事件的概率，或者得到这个随机变随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。量的某些数字特征，并将其作为问题的解。蒙特卡罗方法（统计模拟方法，随机蒙特卡罗方法（统计模拟方法，随机模拟模拟

2、方法）方法）：以概率和统计理论方法以概率和统计理论方法为基础，为基础，使用（伪）随机数来解决计算问题。使用（伪）随机数来解决计算问题。将所将所求解的问题同一定的概率模型相联系，用求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得电子计算机实现统计模拟或抽样，以获得问题的近似解。问题的近似解。模拟通常借助于计算机模拟通常借助于计算机进行：进行：从从理论上讲，理论上讲，蒙特卡罗方法需要大量的实验。计算机技蒙特卡罗方法需要大量的实验。计算机技术的发展，使得蒙特卡罗方法在术的发展，使得蒙特卡罗方法在最近几十最近几十年得到年得到快速的普及。快速的普及。 对于对于那些由于计算过于复杂

3、而难以那些由于计算过于复杂而难以得到解析解或者根本没有解析解的得到解析解或者根本没有解析解的问题，蒙特卡罗方法是一种有效的问题，蒙特卡罗方法是一种有效的求出数值解的方法求出数值解的方法。模拟的作用：模拟的作用： 对对建模过程中的假设进行鉴定，对建模过程中的假设进行鉴定，对理论研究的结论加以检验理论研究的结论加以检验； 对对不同的实现方案进行多次模拟，不同的实现方案进行多次模拟，按照既定的目标函数对不同方案进按照既定的目标函数对不同方案进行比较，从中选择最优方案行比较，从中选择最优方案。设计正确的模拟时间推进机理：模拟过程设计正确的模拟时间推进机理：模拟过程中应根据系统的特性正确推进模拟时间，中

4、应根据系统的特性正确推进模拟时间，使系统中各要素与发生的事件保持同步。使系统中各要素与发生的事件保持同步。通常，模拟时间是模拟的主要自变量通常，模拟时间是模拟的主要自变量。 排队排队系统通常采用离散型系统通常采用离散型模拟模型。模拟模型。其中，发生系统状态变化的事件有两个：其中，发生系统状态变化的事件有两个：一是有顾客到达；二是服务员完成服务。一是有顾客到达；二是服务员完成服务。将最近发生上述两种事件之一的时刻设将最近发生上述两种事件之一的时刻设置为下次事件发生点，就可将服务过程置为下次事件发生点，就可将服务过程描述为下图所示的模拟模型。描述为下图所示的模拟模型。下次下次事件事件法法：应当：应

5、当重点对系统状态可能重点对系统状态可能发生改变的事件进行描述，并确定这些发生改变的事件进行描述，并确定这些事件之间的逻辑关系。事件之间的逻辑关系。初始化数据初始化数据确定顾客到达时间确定顾客到达时间有顾客到达？有顾客到达？服务员空闲服务员空闲确定服务时间确定服务时间是是是是否否否否服务员空闲？服务员空闲？服务员开始服务服务员开始服务顾客进入顾客进入队列等候队列等候推移时间至下次事件发生点推移时间至下次事件发生点队列中有顾客？队列中有顾客？否否是是选择一位顾客选择一位顾客该顾客离开队列该顾客离开队列排队系统中服务过程的模拟模型排队系统中服务过程的模拟模型模拟的一般模拟的一般步骤：步骤：1) 构造

6、构造和描述概率过程和描述概率过程：对于本身就具有：对于本身就具有随机性质的问题，主要是正确描述和模拟随机性质的问题，主要是正确描述和模拟这个概率这个概率过程；对于过程；对于本身不是随机性质的本身不是随机性质的确定性问题，比如计算定积分，必须事先确定性问题，比如计算定积分，必须事先构造一个人为的概率过程，它的某些参量构造一个人为的概率过程，它的某些参量正好是所求问题的解，即将不具有随机性正好是所求问题的解，即将不具有随机性质的问题转化为随机性质的问题质的问题转化为随机性质的问题。最简单、最基本、最重要的概率分布是最简单、最基本、最重要的概率分布是（0，1）上的均匀分布上的均匀分布（或矩形或矩形分

7、布）。随机数就是具有分布）。随机数就是具有这种均匀分布的随机变量，随机数序列就是具有这种均匀分布的随机变量，随机数序列就是具有这种分布的总体的一个简单子样，即一个具有这这种分布的总体的一个简单子样，即一个具有这种分布的相互独立的种分布的相互独立的随机变数序列随机变数序列。产生随机数。产生随机数的问题，就是从这个分布的抽样问题。的问题，就是从这个分布的抽样问题。2) 实现从已知概率分布抽样实现从已知概率分布抽样：构造了概率：构造了概率模型后，由于各种概率模型都可看作是由模型后，由于各种概率模型都可看作是由不同的概率分布构成，因此生成已知概率不同的概率分布构成，因此生成已知概率分布的随机变量，就成

8、为实现蒙特卡罗模分布的随机变量，就成为实现蒙特卡罗模拟（随机抽样实验）的基本手段拟（随机抽样实验）的基本手段。3) 建立各种估计量建立各种估计量：一般说来，构造了：一般说来，构造了概率模型并能从中概率模型并能从中抽样（即实现模拟实抽样（即实现模拟实验）后，验）后，我们就要确定一个随机变量，我们就要确定一个随机变量，作为作为所求所求的问题的的问题的解。即针对模拟实验解。即针对模拟实验的结果考察其统计特性的结果考察其统计特性（样本均值、方（样本均值、方差、置信区间等差、置信区间等），建立），建立各种各种估计量，估计量，从中得到问题的解从中得到问题的解。明确问题，建立明确问题，建立模型模型收集收集和

9、整理数据和整理数据资料资料编制程序编制程序，模拟，模拟运行运行分析分析模拟输出结果模拟输出结果模拟的优点：简单、模拟的优点：简单、快速、适应性强快速、适应性强模拟的缺点模拟的缺点： 建立建立和运行模拟模型可能相当昂贵和运行模拟模型可能相当昂贵； 模拟模型模拟模型的随机性使得结论受到的随机性使得结论受到限制限制能能相对容易地近似很复杂的相对容易地近似很复杂的随机系统，问题随机系统，问题的几何的几何形状的复杂性形状的复杂性对对其其影响影响不大不大；可以可以在广泛的条件下估计候选方案的性能在广泛的条件下估计候选方案的性能；建模建模者可以在不同层次的水平上进行者可以在不同层次的水平上进行控制；控制；产

10、生规则：产生规则：1()mod( ),nnxaxbc特点：特点：序列最终陷入循环序列最终陷入循环( (循环周期至多循环周期至多c c个数个数) )序列各数之间缺乏统计独立性序列各数之间缺乏统计独立性c c一般取得较大，如：一般取得较大，如：15311(2)mod(21).nnxx2 2、随机数、随机数和随机变量的生成和随机变量的生成2.1 2.1 均匀分布随机数的生成均匀分布随机数的生成线性同余法线性同余法(Linear Congruence)(Linear Congruence)2.2 2.2 MatlabMatlab随机模拟函数随机模拟函数1)1)生成生成mxn阶阶0,1均匀分布的随机数矩

11、阵均匀分布的随机数矩阵rand(m,n)( , )normrndm n 211( )exp,.22xf xxR randn(m,n) 产生产生标准正态分布的随机数标准正态分布的随机数2.3 2.3 其它随机变量的生成其它随机变量的生成1）逆变换法）逆变换法1( ),F1()( )( )( ),PxP FxPF xF x 令令逆变换法的具体步骤：逆变换法的具体步骤：确定随机变量确定随机变量X的概率分布函数的概率分布函数F(x)；例：指数分布的分布函数为：例：指数分布的分布函数为：1,0,( )0,0.xexF xx ( ) 1,0 xF xex 1ln(1),x解得解得1ln .x可取可取2）一

12、般离散型分布问题）一般离散型分布问题实际应用中，许多随机变量无法用某种规实际应用中，许多随机变量无法用某种规则的分布来近似，此时我们可以通过累积则的分布来近似，此时我们可以通过累积分布进行处理。分布进行处理。(),0,1,iiPxpi 00,0,()(),0.iiijjiF xPxpi 则则1()(),1,2,iiF xF xi 0,00.3,1, 0.30.5,(0,1).2, 0.50.9,3,0.91,rrXr Urr 模拟方法：模拟方法：例：例：X0123概率概率0.30.20.40.1随机变量函数or 累积分布构造一个构造一个矩形区域包含曲边梯形，矩形区域包含曲边梯形，其其高高h满足

13、满足max( )a x bhf x 由大数定律，事件由大数定律，事件发生的频率发生的频率依概率收敛于事件依概率收敛于事件发生的概率。发生的概率。3、随机、随机模拟案例模拟案例在矩形区域内产生在矩形区域内产生n（足够大）个点（足够大）个点具体思路：具体思路：落在曲边梯形内的点为落在曲边梯形内的点为m个，则个，则则所求定积分则所求定积分( ),baf x dxmpSn( )() .bamf x dxba hn随机点随机点数目数目面积的近面积的近似值似值随机点随机点数目数目面积的近面积的近似值似值1001000.32000.3200100010000.34000.34002002000.39000.

14、390010000100000.33200.33203003000.26670.266750000500000.336760.336765005000.28800.28801000001000000.335860.335865005000.31200.31201000001000000.329560.32956求定积分求定积分 的的Monte-CarloMonte-Carlo模拟结果模拟结果12013x dx 例例2 2：估算圆周率：估算圆周率1-1-11模拟产生在正方形模拟产生在正方形ABCD中均匀分布中均匀分布的的 n个点，如果有个点，如果有m个在圆内，则圆的面积与个在圆内，则圆的面积与正

15、方形面积之比可近似正方形面积之比可近似为为 m/n，即：，即：4.4mmnn模拟程序：估算圆周率模拟程序：估算圆周率n=10000; xy=-1+2*rand(2,n);x=xy(1,:); y=xy(2,:);mypi=4*sum(x.2+y.2=1)/n较简洁实现：较简洁实现：模拟运行结果：模拟运行结果：模拟次数模拟次数n n100003.144810000315201000003.1406例例3 3：求解：求解非线性规划非线性规划min( )( )0,1,2,. .,1,2, .nxRijjjfxgxims taxbjn在估计的区域在估计的区域12(,) |,1

16、,2,njjjx xxxabjn内随机取若干试验点，从试验点中找出内随机取若干试验点，从试验点中找出可行点，再从可行点中选择最小点。可行点，再从可行点中选择最小点。例：例：22121212121212max283310. .10,0zxxx xxxxxstxxxx Matlab中编程需设计三个中编程需设计三个M文件：主程序，文件：主程序，目标函数程序，约束条件测试程序目标函数程序，约束条件测试程序例例4 4：设有一只兔子和一匹狼，兔子在点：设有一只兔子和一匹狼，兔子在点(0,0)(0,0)，狼在点，狼在点(b,0)(b,0)。如果兔子与狼同时。如果兔子与狼同时发现对方，并开始一场追逐。兔子沿发

17、现对方，并开始一场追逐。兔子沿y y轴轴向位于向位于(0,a)(0,a)的巢穴跑，而狼则在其后追的巢穴跑，而狼则在其后追赶（狼追赶的方向应始终正对兔子）。假赶（狼追赶的方向应始终正对兔子）。假设狼的速度是兔子的两倍，且都匀速奔跑，设狼的速度是兔子的两倍，且都匀速奔跑，试确定兔子能否安全跑回巢穴？试确定兔子能否安全跑回巢穴？AB火车运行方向火车运行方向某人某人他已了解到：他已了解到： 1) 火车火车从从A站到站到B站的运行时间是均值为站的运行时间是均值为30分钟，分钟，标准差为标准差为2分钟的随机变量；分钟的随机变量； 思考：思考：请研究他能否赶上这趟火车。请研究他能否赶上这趟火车。例例5：赶上

18、火车的概率赶上火车的概率【问题问题】如图，一列火车从如图，一列火车从A站开往站开往B站，站，某人每天赶往某人每天赶往B站上这趟火车站上这趟火车.出发时刻出发时刻午后午后1:00 午后午后1:05 午后午后1:10频频 率率 0.7 0.2 0.1他到达他到达B 站的时刻的频率分布为站的时刻的频率分布为时刻时刻 午后午后1:28 午后午后1:30 午后午后1:32 午后午后1:34频率频率 0.3 0.4 0.2 0.1他能否及时赶上火车？他能否及时赶上火车？ 2) 火车火车在下午大约在下午大约1 1点离开点离开A站，离开站，离开时刻的时刻的频率频率分布如下：分布如下：明确问题：明确问题：他能及

19、时赶上火车的概率是多少？他能及时赶上火车的概率是多少？i）分析法：）分析法：用概率统计知识用概率统计知识建立模型建立模型，求，求解析解析解解. .ii）模拟法：用概率统计知识建立模型，通过模拟）模拟法：用概率统计知识建立模型，通过模拟求求近似解近似解。 即先建立模拟模型，然后通过计算机模拟得到即先建立模拟模型，然后通过计算机模拟得到问题的近似解。在同样条件下多次试验，计算问题的近似解。在同样条件下多次试验，计算他能及时赶上火车的频率。他能及时赶上火车的频率。建模方向（思路）：建模方向（思路）：问题分析：能及时赶上火车的充要条件是问题分析：能及时赶上火车的充要条件是：其中其中 T1火车从火车从A

20、站出发的时刻；站出发的时刻；T2火车的运行时间；火车的运行时间；T3他到达他到达B站的站的时刻时刻.是什么变是什么变量？如何量？如何模拟？模拟？213213TTTTTT即ii）将午后将午后1时记时记为为 t=0，设火车运行时间，设火车运行时间T2服从正态分布：服从正态分布： T2 N（30，22）.i）假设）假设T1，T2，T3都是相互独立的随机变量都是相互独立的随机变量;基本基本假设：假设：T1（分） 0 5 10P(t) 0.7 0.2 0.1T3（分） 28 30 32 34P(t) 0.3 0.4 0.2 0.1建立模型：为了简化计算，将下午建立模型：为了简化计算，将下午1点记为初始点

21、记为初始时刻。得到随机变量时刻。得到随机变量T1 (火车火车出发出发时刻时刻) 和和 T3 (人人到达到达B站站时刻时刻) 的的分布分布律分别为：律分别为：能及时赶上火车的概率能及时赶上火车的概率 p=PT3T1T2 .如果如果r为在为在(0,1)均匀分布的随机数，为了模拟随机均匀分布的随机数，为了模拟随机变量变量T1和和T3，可以通过如下方法：，可以通过如下方法：. 19 . 0,10, 9 . 07 . 0, 5, 7 . 00, 01rrrt. 0 . 19 . 0,34, 9 . 07 . 0,32, 7 . 03 . 030, 3 . 00,283rrrrt，则则t1和和t3可分别可

22、分别用来模拟随机变用来模拟随机变量量T1和和T3主要变量说明：主要变量说明：n 模拟次数模拟次数k 临时变量，存储当前累计模拟次数临时变量，存储当前累计模拟次数count 存储赶上火车的次数存储赶上火车的次数模拟算法设计模拟算法设计输入：输入：赶火车次数（天数）赶火车次数（天数）输出：输出：赶上火车的频率赶上火车的频率第第 1 步步 输入模拟次数输入模拟次数n第第 2 步步 k=1，count=0第第 3 步步 当当k T3, count=count+1，END第第 10 步步 执行第执行第3步步第第 11 步步 输出赶上火车输出赶上火车频率频率p=count/n算法（分步骤描述）算法（分步骤

23、描述）算法（伪代码描述）算法（伪代码描述）i）初始化：）初始化： 输入模拟次数输入模拟次数n; count=0;ii）模拟）模拟n次次for i=1 to n,模拟随机变量模拟随机变量T1,T2,T3， 分别赋给分别赋给t1,t2,t3;if t1+t2 t3,count=count+1end if end forapp_prb=count / n;n=input(输入模拟次数：输入模拟次数：);count=0;for i=1:n, rt1=rand; %模拟随机变量模拟随机变量t1（火车从（火车从A站出发的时刻）站出发的时刻） if rt1=0.7 & rt10.9 T1=5; el

24、se T1=10; end T2=30+randn*2; %模拟随机变量模拟随机变量t2（火车的运行时间）（火车的运行时间） %模拟随机变量模拟随机变量t3（他到达（他到达B站的时刻）站的时刻） rt3=rand; if rt3=0.3 & rt3=0.7 & rt30.9 T3=32; else T3=34; end if T3 n打印输出打印输出产生装卸服务时间产生装卸服务时间 tutatf统计卸位空闲时间统计卸位空闲时间ti=ti+ta-tf统计船舶等候时间统计船舶等候时间tw=tw+tf-ta设置船舶离港时间设置船舶离港时间tf=ta+tu设置船舶离港时间设置船舶离港时

25、间tf=tf+tu是是是是否否否否k=k+1模拟模拟次数次数到达船到达船舶数舶数模拟总时间模拟总时间（小时）（小时）船舶平均船舶平均滞留时间滞留时间（小时）（小时）船舶滞留船舶滞留最大时间最大时间（小时）（小时）滞留队列滞留队列最大长度最大长度装卸泊位装卸泊位利用率利用率110007860752.42300586%210008106737.85274483%310008039345.80319584%410007950642.30247484%510007864543.92289586%平均平均10007964444.462864.684.6%表：港口对船舶服务的模拟结果表：港口对船舶服务的模

26、拟结果到达间隔时间：到达间隔时间：15,14515,145，装卸时间：，装卸时间：45,9045,90模拟模拟次数次数到达船到达船舶数舶数模拟总时间模拟总时间（小时）（小时）船舶平均船舶平均滞留时间滞留时间（小时）（小时）船舶滞留船舶滞留最大时间最大时间（小时）（小时）滞留队列滞留队列最大长度最大长度装卸泊位装卸泊位利用率利用率110008092814.33173368%210007904414.84144370%310007894816.42167370%410008206111.61138367%510008055814.60123370%平均平均10008030814.36149369%

27、表：港口对船舶服务的模拟结果表：港口对船舶服务的模拟结果到达间隔时间：到达间隔时间：15,14515,145，装卸时间：，装卸时间：35,7535,75模拟模拟次数次数到达船到达船舶数舶数模拟总时间模拟总时间（小时）（小时）船舶平均船舶平均滞留时间滞留时间（小时）（小时）船舶滞留船舶滞留最大时间最大时间（小时）（小时）滞留队列滞留队列最大长度最大长度装卸泊位装卸泊位利用率利用率110006461739.09274585%210006431439.62380785%310006514648.28307685%410006409846.89247586%510006599331.12226582%

28、平均平均10006483441.002875.684.6%到达间隔时间：到达间隔时间：10,12010,120，装卸时间：，装卸时间：35,7535,75表：港口对船舶服务的模拟结果表：港口对船舶服务的模拟结果到达间隔到达间隔( (小时小时) )56 7 89101112131415161718到达船舶到达船舶数数( (艘艘) )13 6 791011111197654对进入该港口的对进入该港口的100100艘船进行实际调查，艘船进行实际调查，记录其活动情况，得到这记录其活动情况，得到这100100艘船到达港艘船到达港口的时间间隔和装卸时间的分布口的时间间隔和装卸时间的分布情况情况. .装卸时

29、间（小时）装卸时间（小时）910111213141516装卸船舶数（艘）装卸船舶数（艘）2022191610832表：表：100100艘船到达港口的时间间隔频数表艘船到达港口的时间间隔频数表表：表：100100艘船装卸时间频数表艘船装卸时间频数表到达间隔56789101112131415161718累积频率0.010.040.100.170.260.360.470.580.690.780.850.910.961.00装卸时间910111213141516累积频率0.200.420.610.770.870.950.981.00表：船舶到达港口的时间间隔累积分布表表：船舶到达港口的时间间隔累积分布

30、表表：船舶装卸时间累积分布表表：船舶装卸时间累积分布表为为比较准确地反映系统的性能，我们应采用比较准确地反映系统的性能，我们应采用稳态模拟方式，可以考察系统稳态模拟方式，可以考察系统运行运行 360360天天的的情况。假定港口每天情况。假定港口每天2424小时连续工作，模拟小时连续工作，模拟时间时间T T设置设置为为 86408640小时小时. .考虑考虑到船舶到达港口的事件是影响整个系统到船舶到达港口的事件是影响整个系统状态的主要因素，可以将模拟时间设置为该状态的主要因素，可以将模拟时间设置为该事件的发生时刻事件的发生时刻。初始化数据初始化数据 t=tf=tw=ti=0产生船舶到达间隔时间产

31、生船舶到达间隔时间 tb设置模拟时间设置模拟时间 t=t+tbt8640打印输出打印输出产生装卸服务时间产生装卸服务时间tut= 360?是是否否w=w+1统计前统计前w个班的数据个班的数据w=20？是是否否卸位空闲？卸位空闲？否否是是统计卸位空闲时间统计卸位空闲时间统计等待卸车时间统计等待卸车时间设置开始卸车时间设置开始卸车时间计算卸完时刻计算卸完时刻Cc.t, tbuCc.s=4卡车数卡车数（辆）（辆）每班电铲平每班电铲平均空闲时间均空闲时间（分钟）（分钟）每班每辆卡车每班每辆卡车装车等待时间装车等待时间（分钟）（分钟）每班每辆卡车每班每辆卡车卸车等待时间卸车等待时间（分钟）（分钟）电铲效

32、率电铲效率每辆卡车每辆卡车平均效率平均效率平均班产平均班产量（吨）量（吨）5139.97.611.861.14%94.62%3633.46102.010.817.271.66%92.23%4253.6771.115.125.780.25%88.67%4757.9847.822.236.586.73%83.69%5120.5934.730.949.390.36%77.74%5334.21027.238.463.692.45%71.67%5453.11125.342.980.392.98%65.78%5480.8表：不同卡车数的模拟结果表：不同卡车数的模拟结果图：卡车数量与电铲和卡车利用率的关系图

33、：卡车数量与电铲和卡车利用率的关系可以看出，电铲效率随卡车数可以看出，电铲效率随卡车数m的增加而的增加而上升，而卡车效率却随上升，而卡车效率却随m的增加而减少。的增加而减少。若要求电铲的效率在若要求电铲的效率在85%以上，则应保证以上，则应保证m7。此时，为保证卡车有较高效率。此时，为保证卡车有较高效率（75%），可取），可取m=8或或9.从从班产量的角度看，电铲效率越高（即班产量的角度看，电铲效率越高（即卡车数卡车数m越大），班产量相应越高，因此越大），班产量相应越高，因此应取应取m7。从表中看出，将卡车数从。从表中看出，将卡车数从9辆辆提高到提高到10辆时，班产量只提高了约辆时，班产量只提

34、高了约2.23%，却付出了卡车效率降低却付出了卡车效率降低6.07%的代价，所的代价，所以最多只考虑以最多只考虑9辆卡车的辆卡车的情况情况.4、模拟、模拟结果的统计分析结果的统计分析 在在对一个具有随机因素的系统设计出计算对一个具有随机因素的系统设计出计算机模拟模型后，不能只将单次模拟运行结果机模拟模型后，不能只将单次模拟运行结果看作问题的真实答案。模拟的输出结果是分看作问题的真实答案。模拟的输出结果是分布特征未知的随机变量，每次运行的结果仅布特征未知的随机变量，每次运行的结果仅仅是对该随机变量所有观察值总体的一次抽仅是对该随机变量所有观察值总体的一次抽样，对总体的代表性不够，虽然可以通过增样，对总体的代表性不够，虽然可以通过增加模拟运行的时间从而增加抽样次数，但这加模拟运行的时间从而增加抽样次数，但这些数据总是由一个些数据总是由一个“种子种子”经过一定的算法经过一定的算法而获得的伪随机序列，它们是自相关的，并而获得的伪随机序列，它们是自相关的，并不能构成统计上独立的随机样本。不能构成统计上独立的随机样本。1）系统性能测度）系统性能测度考察系统的性能，就是根据随机样本的考察系统的性能，就是根据随机样本的取值来估计系统真实参数的取值来估计系统真实参数的统计计量统计计量。一般可用点估计和区间估计。一般