《数学实验》第16讲
主要内容：
随机现象及模拟
随机变量模拟(随机模拟的基础) 蒙特卡罗方法原理
应用实例
随机现象及模拟
一、确定性现象
在一定条件下，某种结果必然会发生的现象
比萨 斜塔 试验
在标准大气压下，水 加热到100摄氏度，
就必然会沸腾。
确定性现象的特点：
可事前预言
二、不确定性现象
小球将落
硬币 将出 现哪 一面?
随机现象及模拟
不确定性现象的特点：不可事前预言
但在大量重复试验，某些不确定现象又呈现出规律性。
随机现象及模拟
称大量同类随机现象所呈现的固有规律为随机现象的统 计规律性。
三、随机现象
为了探索随机现象的规律性，理想化的方法是在相同条 件下将实验大量重复进行，采集到试验数据，再对数据 进行统计分析，得到其规律性。
但当试验周期长，或一个试验具有破坏性时，通过试验 采集数据是不可能的，此时，通过计算机做随机模拟的
随机现象及模拟 四、随机模拟
小球将落入哪 一格?
引例 高尔顿钉板试验
随机现象及模拟
需建立数学模型，描述小球的运动过程
2.小球入口处水平位置为坐标原点0;
假设：
4.小球在不同层向左或向右是相互独立的．
3.小球在每层碰到钉子后，向左或向右等可能位移 一格，不会出现跳格（位移2格以上）的情况；
1.共有n层钉子;
随机现象及模拟
1,
k 1, X
= −
第k层向右 第k层向左
完整描述了小球在n层钉板的运动过程．
k= 1,…,n
则随机变量序列：－1 1 X
kp 1/2 1/2
1 n
n k
k
Y X
=
=
关注：小球经n次碰撞后在钉板底层所处位置
随机变量模拟(随机模拟的基础)
回顾：
高尔顿钉板试验中，小球最终的位置
1 n
n k
k
Y X
=
=
其中
X
k－1 1
p 1/2 1/2
要模拟小球的运动轨迹，首先要模拟随机变量Xk，那么如何模拟 随机变量呢？
一、随机数的生成
函数名 解释
rand 生成(0,1)区间上均匀分布的随机数 unifrnd 生成指定区间内均匀分布的随机数
randn 生成服从标准正态分布的随机数
normrnd 生成指定均值、标准差的正态分布随机数 exprnd 生成服从指数分布的随机数
注：上述命令常被称为伪随机数生成器。
随机变量模拟(随机模拟的基础)
rand(m) 生成m*m维的随机数 rand(m,n) 生成m*n维的随机数
rand([m,n,p ...]) 生成排列成m*n*p... 多维向量的随机数
基本语法如下：
问题: 如何模拟在区间[a, b]内均匀分布随机数？
1. a+(b-a)*rand(m, n) 2. unifrnd(a, b, m, n)
随机变量模拟(随机模拟的基础)
二、离散型随机变量
思考：如何利用rand生成1000个下列离散型随机变量? （5min练习）
分析：
rand是生成（0,1）上均匀分布随机数，生成数落在 (0,0.5)和[0.5,1)上概率均为0.5，故可令－1 1 X
kp 1/2 1/2
=
5 . 0 1
5 . 0 1
rand ,
rand ,
X
k-
随机变量模拟(随机模拟的基础)
参考程序：
N=1000;
Y=rand(1,N);
for
i=1:N
if
Y(i) < 0.5 X(i)= -1;
else
X(i)= 1;
end end
X
思考：一般的离散随 机变量如何模拟？
随机变量模拟(随机模拟的基础)
蒙特卡罗方法原理
引例: 请计算
分析 dx 原函数不存在，所以无法用Newton-Leibniz公式求解。
x (x) I =
sindx f
b
a(x)
x b a n
= −
x
i= + a i x i ( = 0,1, , ) n 故可考虑数值积分方法：利用几何意义计算
将分成 [ , ] a b n等分:
x dx (x) I
2sin
1=
记 f x ( )
i= y
i( i = 0,1, , ) n ( ) d
b
a
f x x
+ + + y x y x
0 1y
n−1 x
其他近似计算方法：
思路：定积分几何意义为曲边梯形的面积，可考虑求曲边梯形与矩形
0 1 1
(
n)
b a y y y
n
−= − + + +
(左矩形公式) 蒙特卡罗方法原理
问题：如何求面积之比？
可考虑把面积之比看成概率，于是可用频率去近似。因此可在矩阵 区域内均匀投点，求解两个区域内点数之比。
基本思想：
当所要求解的问题是某种事件出现的概率，或者是某个随机变量 的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现 的频率，或者这个随机变数的平均值，并用它们作为问题的解
蒙特卡罗方法原理
蒙特卡罗方法
蒙特卡罗方法，或称计算机随机模拟方法，是一种基于“随机 数”的计算方法。源于美国在第二次世界大战研制原子弹的
“曼哈顿计划”，该计划的主持人之一数学家冯·诺伊曼用驰名 世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙 上了一层神秘色彩。
蒙特卡罗方法原理
分析：
设圆的半径为1，则其外切正方形的面 积为4，圆的面积为 ，现在模拟产生 在正方形ABCD中均匀分布的点n个，
如果这n个点中有m个点在该圆内，则 圆的面积与正方形ABCD的面积之比 可近似为m/n；即
例： 采用蒙特卡罗方法估计圆周率
n m n
m
44
思考如何编程？(5min) 蒙特卡罗方法原理
程序参考代码：
n=input('请输入产生点的个数：');
m=0;
for i=1:n
x=-1+2*rand; y=-1+2*rand;
if x^2+y^2<=1 % if within the circle m=m+1; % count 1
end end
mypi=4*m/n
蒙特卡罗方法原理
应用实例
例1: 曲线所围区域面积 （5min）
用随机模拟方法估算两条抛物线y =x2 , x = y2所围图形的面积
分析:
(1)两条曲线交点, [0,0], [1,1]
(2)用一个矩形区域包含所围 区域:
(3)所围区域内的点满足:
x
2 y x
1 0
, 1
0 x y
应用实例
N=1000;
data=rand(N,2);
x=data(:,1);
y=data(:,2);
II=find(y<=sqrt(x)&y>=x.ˆ2);
M=length(II);
S=M/N
曲线所围区域面积---程序
运行输出：
S =
0.3276
应用实例
例2: 求解下列优化模型
− +
−
= 5
. 3 .
) 2 (
4 )
1 (
3 )
, ( min
2 1
2 2
2 1
2 1
x t x s
x x
x x f
分析：可在可行域内随机投点，比较各点处函数值的大小找到最小 值
function
y=fun(x)y=3*(x(1)-1).^2 + 4*(x(2)-2).^2;
应用实例
function opt_sim_ex1 n = 1e5; %随机点个数 fobj = inf;
for i = 1:n,
x(1) = 3*rand; x(2) = 5*rand;
temp = fun(x);
if temp < fobj, %找"最小值"
fx = x; fobj = temp;
end end
%返回近似最优解
fx %变量x1,x2取值 fobj %对应函数值
更多例子（回到第1讲）：
限定区域的随机投点实验
应用实例（补充）
(一) 二重积分几何意义是计算体积: 由于D的边界曲线交点为
：(1,-1), (4,2), 被积函数在求积区域内的最大值为16。积分值 是一个三维图形所围体积(V1) ，该三维图形位于立方体区域 其中D为 y= x – 2 与 y2 = x 所围区域。
计算二重积分：
D
dxdy xy2
例1： 估算二重积分
分析:
(二) 蒙特卡罗方法求解: 向立方体内投N个随机点，统计落在曲 顶柱体内的个数M，则
{(x，y，z) | 0≤ x ≤4，–1≤ y ≤2，0 ≤ z ≤16 }
该立方体区域的体积(V2)为192 应用实例（补充）
function testmain N=100000;
for k=1:7 %多次模拟 V1(k)=mysim(N);
end V1
function v =mysim(N) V2=192; %V2=4*3*16
d =rand(N,3); x =4*d(:,1);
y =-1+3*d(:,2); z =16*d(:,3);
%下面表达式可结合find和length完成
M =sum((x>=y.^2) & (x<=y+2) & (z<=x.*y.^2) );
估算二重积分 —程序 理论结果： 7.5857
V1 =
7.5859 7.5898 7.6262 7.5629 7.7894 7.6781 7.4304
应用实例（补充）
例3： 相遇问题
甲乙两船在24小时内独立地随机到达码头. 设两船到达码头时刻分 别为 X 和 Y，且X ~ U(0 , 24), Y ~ U(0 , 24)。如果甲船到达码头后停 留2小时,乙船到达码头后停留1小时.问两船相遇的概率有多大？
1
且 +
x x y y
相遇条件：
(1)甲比乙先到码头：
(2)乙比甲先到码头：
2
且 +
y y x
x
应用实例（补充）
概率值：P = 0.1207
相遇问题程序
N=2000;
P=24*rand(2,N);
X=P(1,:);Y=P(2,:);
I=find(X<=Y&Y<=X+2);
J=find(Y<=X&X<=Y+1);
运行程序：
F =0.1175
应用实例（补充）
学到了什么？
随机现象及模拟
随机变量模拟(随机模拟的基础) 蒙特卡罗方法原理
应用实例
