第四章 随机随机模拟方法1

①借助区间(0,1)分布产生的随机数，对 变量 T1、 T3概率分布进行统计模拟；
②根据变量 T1、T2 、T3 概率分布及模拟 程序、命令产生n 个随机分布数； ③使用随机产生的n 组随机数验证模型中 的关系表达式是否成立； ④计算n 次模拟实验中，使得关系表达 式成立的次数k ； ⑤当n 时，以

（三）puffon随机投针实验


1777年Puffon（法）提出用投针实验求圆 周率Pi的问题。 间距为a的平行线，随机投掷一枚长为l（ l<a）的针，试求此针与一平行线相交的概 率P。
4
5百度文库
6

（四）Puffon投针的R实现
>source("buffon.R") >buffon(100000,l=0.8,a=1)##直接调用刚 刚编写好的“buffon”函数， 当然，l和a的取值由于在函数编写中已经 给出了，因此这里不必再给出， 可以直接用buffon（100000）。
进行计算机统计模拟的基础是抽象现实系统的 数学模型
为了便于建模，对模型中使用的变量作出如 下假定：
T1 ：火车从A站出发的时刻；
T2 ：火车从A站到B站的运行时间；
：某人到达B站的时刻； ：随机变量 T2 服从正态分布的均值； ：随机变量 T2 服从正态分布的标准差；
T3


T 为了分析简化，假定13时为时刻t=0，则变量 T1 、 3 的分布律为：
12
第三节 用蒙特卡洛方法求解确 定性问题

（一）计算定积分


为了简化计算，a=0，b=1。 计算定积分值也就是求曲边梯形的面积S， 常用方法有： (1).随机投点法
13
14

记录实验次数N，成功次数M，用M/N作为概 率p的估计值，即可得出定积分I的近似解。
15
16

(2).平均估值法
7
另一种求pi的方法
8
9
MC1 <- function(n){ k <- 0; x <- runif(n); y <- runif(n) for (i in 1:n){ if (x[i]^2+y[i]^2 < 1) k <- k+1##在n个值当中，有多少个i满 足if函数 } 4*k/n }
k 作为此人能赶上 n
是
产生随机数 验证模型
关系式 成立 否
成立次数k=k+1
成立次数不变
试验次数 是否达到n次 是 计算估计结果 k/n
否
火车的概率p 的近似估计；
R 程序：
windows(7, 2)#作图窗口大小 prb = replicate(10, { #括号内程序重 复100次 x = sample(c(0, 5, 10), 1, prob = c(0.7, 0.2, 0.1)) y = sample(c(28, 30, 32, 34), 1, prob = c(0.3, 0.4, 0.2, 0.1)) plot(0:40, rep(1, 41), type = "n", xlab = "time", ylab = "", axes = FALSE) axis(1, 0:40) r = rnorm(1, 30, 2) points(x, 1, pch = 15) i=0 while (i <= r) { i=i+1 segments(x, 1, x + i, 1) if (x + i >= y) points(y, 1, pch = 19) Sys.sleep(0.1) } points(y, 1, pch = 19) title(ifelse(x + r <= y, "poor... missed the train!", "Bingo!catched the train!")) Sys.sleep(4) x+r>y }) mean(prb)
>source("MC1.R");MC1(100000) [1] 3.14268
10

（五）统计模拟的一般步骤
11
第二节 随机模拟方法的特点



（一）方法新颖、应用面广、适用性强 （二）随机模拟方法的算法简单，但计算 量大 （三）模拟结果具有随机性，且精度较低 （四）模拟结果的收敛过程服从概率规律 性
第四章 随机模拟方法




第一节 第二节 第三节 第四节 应用 第五节 拟 第六节
概述 随机模拟方法的特点 用蒙特卡罗方法求解确定性问题 随机模拟方法在随机服务系统中的 集装箱专用码头装卸系统的随机模
随机模拟方法在理论研究中的应用
1
第一节 概述

（一）随机（统计）模拟的定义
随机模拟即是计算机统计模拟，它实质上是计算机 建模，而这里的计算机模型就是计算机方法、统计模 型(如程序、流程图、算法等)，它是架于计算机理论 和实际问题之间的桥梁。它与统计建模的关系如下图。
17
18
例：赶火车问题
一列列车从A站开往B站，某人每天赶往B站上车。他 已经了解到火车从A站到B站的运行时间是服从均值为 30min，标准差为2min的正态随机变量。火车大约下午13： 00离开A站，此人大约13:30到达B站。火车离开A站的时刻 及概率如表1所示，此人到达B站的时刻及概率如表2所示。 问此人能赶上火车的概率有多大？
T1 / min
P(t )
0
5
10
T3 / min 28 30 32 34
P(t )
0.7 0.2 0.1
0.3 0.4 0.2 0.1
T 此人能及时赶上火车的充分必要条件为：1 T2 T3 ，所 p 以此人能赶上火车的概率模型为： {T1 T2 T3}。
R软件求解的总算法：
编写R程序
实际问题 统计、逻辑 模型
计算机模拟（程序、算法）
实际解 统计、计算机解
2

（二）随机模拟方法
一般地，随机模拟分类如下： 若按状态变量的变化性质分为连续随机模拟和离散随 机模拟。 而按变量是否随时间变化又可分为动态随机模拟和静 态随机模拟。 常用的随机模拟方法主要有以下几种： 1.蒙特卡罗法 2.系统模拟方法 3.其它方法：包括Bootstrap(自助法)、MCMC（马氏 链蒙特卡罗法）等。
[1]3.142986 还可以更改n的数值，与l、a的数值（见 表4-1） > buffon(100000,l=0.8,a=1) [1] 3.119334 > buffon(1000000,l=0.8,a=1) [1] 3.136941 调用已编号的R程序buffon.R进行模拟， > buffon(10000000,l=0.8,a=1) 取n=100000， [1] 3.142697 l=0.8，a=1。
#应用R软件对buffon投针实验进行模拟 pi的取值 #先编写“buffon”函数， buffon <- function(n, l=0.8, a=1){ k<-0 theta <- runif(n, 0, pi); x <- runif(n, 0, 1/2) for (i in 1:n){ if (x[i]<= l/2*sin(theta[i])) k <- k+1 } 2*l*n/(k*a)##函数最终输出的值，比如 当给定参数n时，输出2*l*n/(k*a) }
24
本节课结束
你学到了什么？
25
表1：火车离开A站的时刻及概率
火车离站时刻 概率
13:00 0.7
13:05 0.2
13:10 0.1
表2：某人到达B站的时刻及概率
人到站时刻 概率
13:28 0.3
13:30 0.4
13:32 0.2
13:34 0.1
——问题的分析——
这个问题用概率论的方法求解十分困难，它 涉及此人到达时刻、火车离开站的时刻、火车 运行时间几个随机变量，而且火车运行时间是 服从正态分布的随机变量，没有有效的解析方 法来进行概率计算。在这种情况下可以用计算 机模拟的方法来解决。