复杂的随机模拟案例

• • • • • • • • • • • • •
If isequal(Freq,[4;4]) Freq0=Freq0+[1 0 0 0 0]; Elseif isequal(Freq,[3;5]) Freq0=Freq0+[0 1 0 0 0]; Elseif isequal(Freq,[2;6]) Freq0=Freq0+[0 0 1 0 0]; Elseif isequal(Freq,[1;7]) Freq0=Freq0+[0 0 0 1 0 ]; Else Freq0=Freq0+[0 0 0 0 1]; End End Em= Freq0*[-3, 0.2,0.5,1,10]’/n;
• 由以上结果可以看到，随着随机抽样次数 的增大，所求概率的模拟值逐渐趋近于理 论值2/3。
CASE2. 抽球问题
• 一袋子中有n个球，从中有放回地抽取 m（n≤m）次，求袋子中的每个球都能 被抽到的概率。这个问题也可以描述为 m（n≤m）个球随机地落到n个盒子中， 求每个盒子中都有球的概率。
数学理论求解
奖金（罚 金）/元
A 8:0 10
B 7:1 1
C 6:2 0.5
D 5:3 0.2
E 4:4 -3
• 注：表中-3表示受罚3元
此游戏从表面上看非常有吸引力，5中可 能出现的结果有4中可得到奖金，且最高 奖金达10元，而只有一种情况受罚，罚 金只有3元。试分析此游戏中谁是真正的 赢家？
数学理论求解
随机模拟方法求解
• 给16个球分别从1-16进行编号，假设8个红 球的编号1-8，而8个白球的编号为9-16。 • 在MATLAB中进行n次随机模拟，每次模拟 用randsample函数生成8个随机整数（取 值范围从1-16）作为一次抽取的8个球。 • 统计n次模拟中各种可能的结果（见下表）出 现的频数 nA , nB , nC , nD , nE ,从而可得摸球者在一次 游戏中得到的奖金（罚金）的数学期望的模拟 值为
P( A1 A2
An )Байду номын сангаас 0
1-
=
袋子中的每个球都能被抽到的概率为
P ( A ) P ( A1 ) P ( A2 ) … P ( An ) 1 P ( A1 A 2 … A n )
1 P( Ai ) ( 1)
i 1 n 2 1 1i j n
• num = 0; • x = 0; % 计算模拟概率 • for i = 1:N • x = randsample(n,m,’true’); % 有放回随机 抽样如果n个数都被抽到，将计数器的值增加1 • if numel(unique(x)) == n • num = num + 1; • end • end • p = num/N; % 模拟概率
MATLAB程序代码如下：
• function p = SheepAndCar(n)
% p = SheepAndCar(n), 利用蒙特卡洛方法求解蒙提霍尔问题，求参赛 者更换选择之后 % 赢得汽车的概率p。这里的n是正整数标量或向量，表示随机抽样的次数。 •
• for i = 1:length(n) • x = randsample(3,n(i),’true’); %随机抽样 • p(i) = sum(x~=3)/n(i); %概率的模拟值 • end
• 设A=“袋子中的每个球都能被抽到”, • Ai =“第i个球没有被抽到”i=1,2,…,n 则有 n 1 m ) • P( Ai ) ( i=1,2,…,n n 2 n )m • P( Ai Aj ) ( i≠j,i,j=1,2,…,n n • P( Ai Aj Ak ) ( n 2 )m i≠j≠k,i,j,k=1,2,…,n n • …… •
• for i = 1:n • x = randsample(16,8,’flase’); % 不放回随机抽样 x(x <= 0) = 1; % 将x中取值为1至8的元素改为1 用来标记红球 x(x>8) = 2; % 将x中取值为9至16的元素改为2用来标记白球 %统计x中1和2出现的次数，整理成[4 4],[3 5],[2 6],[1 7],[0 8] 的形式 • x = sort(x); % 将x从小到大排序 • x1 = diff(x); % 求差分 • x1(end+1)=1; • x1=find(x1); • x1=[0;x1]; • Freq=sort(diff(x1)); • If Freq==8 • Freq=[0;8]; • End
数学理论求解
• 由于游戏开始是参赛者是从三扇门中随机地选 取一扇门，所以在更换选择之前，参赛者赢得 汽车的概率为1/3。 • 经分析可知，若参赛者一开始选中汽车，则更 换选择后一定选不到汽车； • 若参赛者一开始没有选中汽车，则更换选择后 一定能选到汽车。
• 为了求解参赛者更换选择之后赢得汽车的概率， 这里引入两个随机事件： • A= 一开始选中汽车 B=更换选择后选中汽车 根据全概率公式可求得参赛者更换选择之后赢 得汽车的概率为
MATLAB函数probmont代码如下：
% [p0,p] = probmont(n,m,N)，有n个球，从中有放回地抽取m次，求每个 球都能被取到的理论概率p0和蒙特卡洛模拟概率。输入参数N为随机模 拟次数。当抽球次数m小于球的总数n时，理论概率和模拟概率均为0 • • • • • • • • function[p0,p] = probmont(n,m,N) if n>m p0 = 0; % 理论概率 p = 0; % 模拟概率 return; end i = 0:n; % 定义一个向量计算理论概率 p0 = sum((-1).^i*factorial(n)./(factorial(i).*factorial(ni)).*(1-i/n).^(m));
CASE3. 街头骗局揭秘
• 街头常见一类“摸球游戏”，游戏规则： • 一袋中装有16个大小、形状相同的玻璃 球，其中8个红色、8个白色。 • 游戏者从中一次摸出8个，8个球中，当 两种颜色出现以下比数时，摸球者可得到 相应的“奖励“或”惩罚“，如表所列
表：摸球游戏的5中可能结果及奖金（罚金） 可能的 结果
• 设摸球者在一次游戏中得到的奖金（罚金）为X， 则 P( X 10) P( A) 2 0.0001554
1 2C8 C8 P( X 1) P( B ) 0.009946 6 C18 6 2 2C8 C8 P( X 0.5) P(C ) 0.1218 6 C185 3 2C8 C8 P( X 0.2) P( D ) 0.4873 6 C18
• 假设n=20,m=50,也就是说有20个球，从中有放回地随机抽 取50次，则可调用probmont函数计算每个球都能被抽到的 理论概率和模拟概率。 • >> [p0,p] = probmont(20,50,10) % 模拟10次 • p0 = 0.1642 % 理论概率 • p = 0.2000 % 模拟概率 • >>[p0,p] = probmont(20,50,1000) %模拟1000次 • p0 =0.1642 • p =0.1690 • >>[p0,p] = probmont(20,50,10000) % 模拟10000次 • p0 =0.1642 • p =0.1640
这个游戏的玩法是：
参赛者面前有三扇关闭的门，其中一扇门的后 面藏有一辆汽车，而另外两扇门的后面则各藏 有一只山羊。参赛者从三扇门中随机选取一扇， 若选中后面有车的那扇门就可以赢得该汽车。 当参赛者选定了一扇门，但尚未开启它的时候， 节目主持人会从剩下两扇门中打开一扇藏有山 羊的门，然后问参赛者要不要更换自己的选择， 选取另一扇仍然关上的门。这个游戏涉及到的 问题是：参赛者更换自己的选择是否会增加赢 得汽车的概率？
1 2 2 P( B ) P( A) P( B | A) P( A) P( B | A) 0 1 3 3 3
参赛者更换选择后赢得汽车的概率增大了，从最 初的1/3变为2/3了，显然参赛者应该更换自己的 选择。
随机模拟方法求解
• 设两只羊的编号分别为1和2，汽车的编号为3。 • 现在从数字1、2、3中随机选取一个数字，若一 开始选中1或2，则更换选择后选中3，即赢得汽 车；若一开始选中3，则更换选择后选中1或2， 即得不到汽车。 • 将这样的试验重复进行n次，记录一开始选中1 或2的次数m（即更换选择后赢得汽车的次数）， 从而可以确定更换选择后赢得汽车的频率m/n。 由大数定律可知当试验次数n增大时，频率m/n 趋近于更换选择后赢得汽车的概率。
6 C18 7
C84C84 P( X 3) P( E ) 6 0.3807 C18
可得X的分布列
X P 10 1 0.5 0.2
0.4873
-3
0.3807
0.0001554 0.009946 0.1218
• 所以可得X的数学期望（均值）为 E(X) = 10 P(X = 10) + 1 P(X = 1) + 0.5 P(X = 0.5) + 0.2 P(X = 0.2) - 3 P(X = -3) = -0.9732(元) • 即从平均意义上来说摸球这在一次游戏中 要赔掉0.9732元，故此游戏中庄家（游 戏经营者）才是真正的赢家。
10n A nB 0.5nC 0.2nD 3nE Em n
MATLAB函数GamrMont1代码
• function [Em,E0] = GameMont1(n)% 摸球游戏揭秘 • % [Em,E0] = GameMont1(n),求摸球者在一次摸球游戏中得 到的奖金（罚金）的数学期望（均值）的理论值E0和蒙特卡洛模 拟值Em。输入参数n为游戏的次数，它是一个整 数。 • a = nchoosek(16,8); % 组合数（16选8） • p = 0; % 通过循环计算每种可能结果的理论概率 • for i = 4:8 • p(i – 3) = 2^(i~=4) * nchoosek(8,i) * nchoosek(8,8-i)/a; • end • E0 = p * [-3,0.2,0.5,1,10]’; % 数学期望的理论值 • Freq0 = zeros(1,5); % 做n次摸球游戏，计算各种可能的结 果出现的频数 •
第一部分 模拟与概率
肖柳青 主讲 lucyxiao@sjtu.edu.cn
PUB:SSMA_xiao@yeah.net
第六章 复杂的模拟举例
[1] 有趣的蒙提霍尔问题（Monty Hall problem） [2] 抽球问题 [3] 街头骗局揭秘 [4] 求圆周率π(另一种用蒙特卡洛方法) [5] 四人追逐问题 [6] 排队系统模拟实例
CASE1. 有趣的蒙提霍尔问题 （Monty Hall problem）
• 蒙提霍尔问题（Monty Hall problem），也 称为三门问题，
• 是一个源自博弈论的数 学游戏问题，问题的名 字来自美国的电视游戏 节目：Let’s Make a Deal，该节目的主持人 名叫蒙提· 霍尔（Monty Hall）
• SheepAndCar 函数代码的注释部分给 出了该函数的调用格式。
• 下面调用SheepAndCar函数，针对不同的n，求参 赛者更换选择之后赢得汽车的概率的模拟值 • >> p = SheepAndCar([10,100,1000,10000,100000]) %求概率模拟值向量 • p= • 0.7000 0.6600 0.6650 0.6600 0.6663 0.6666

n
P( Ai Aj )
m n
( 1) n 1 P( A1 A2
m
An )
n i i i n i 1 ( 1) C ( 1) Cn n n i 1 i 0
n i 1 i n
随机模拟方法求解
• 首先给n个球从1~n分别编号。 • 在MATLAB中进行N次随机模拟，每次模拟用 randsample函数（或randint函数）生成m 个随机整数（取值范围从1~n）作为m次有放 回抽球，如果这m个随机整数包含了全部的n 个编号，则将计数器的值加1，这样就可以计 算出N次模拟中n个球都能被取到的概率。 • 随着N的增大，这个频率就会越来越接近于每 个球都能被抽到的概率P(A)。