随机行为的模拟:随机抛掷硬币和骰子出现特定面的概率——蒙特卡罗方法的计算机模拟1摘要对蒙特卡罗(Monte Carlo)方法的简介并概述了蒙特卡罗方法的概念、应用领域、求解步骤。
以抛掷硬币和骰子为例,论述了蒙特卡罗方法模拟随机行为的基本思想和基本原理。
给出了实现计算机模拟的MATLAB程序,并且通过最高达千万次级别的计算机模拟试验,准确地模拟了随机抛掷硬币和骰子出现特定面的概率。
2关键词蒙特卡罗(Monte Carlo)方法方法;计算机模拟;随机行为;模拟;概率;MATLAB 程序3引言3.1蒙特卡罗方法的概述:蒙特卡罗方法又称统计模拟法、随机抽样技术,是一种随机模拟方法,以概率和统计理论方法为基础的一种计算方法,是使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。
将所求解的问题同一定的概率模型相联系,用电子计算机实现统计模拟或抽样,以获得问题的近似解。
为象征性地表明这一方法的概率统计特征,故借用赌城蒙特卡罗命名。
3.2蒙特卡洛模拟法简介:蒙特卡洛(Monte Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。
具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时,可用随机模拟法近似计算出系统可靠性的预计值;随着模拟次数的增多,其预计精度也逐渐增高。
由于涉及到时间序列的反复生成,蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的,因此只是在近些年才得到广泛推广。
3.3 蒙特卡洛模拟法提出:蒙特卡罗方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。
数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo —来命名这种方法,为它蒙上了一层神秘色彩。
在这之前,蒙特卡罗方法就已经存在。
1777年,法国Buffon 提出用投针实验的方法求圆周率。
3.4 蒙特卡洛模拟法的应用领域:(1)、直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统,得到某些参数或重要指标。
(2)、蒙特卡洛积分:利用随机数列计算积分,维数越高,积分效率越高。
(3)、MCMC:这是直接应用蒙特卡洛模拟方法的推广,该方法中随机数的产生是采用的马尔科夫链形式。
(4)、蒙特卡罗方法在金融工程学,宏观经济学,生物医学,计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。
3.5 蒙特卡罗解题归结为三个主要步骤:(1)、构造或描述概率过程;(2)、实现从已知概率分布抽样; (3)、建立各种估计量。
4 问题重述蒙特卡罗模拟的真正威力在于对随机行为建模。
从长期来看,一个事件的概率可以视为比值:事件的总数有效的事件数概率 )(A P下面3个随机模型:(1)、抛掷一枚正规的硬币 (2)、抛掷一个正规的骰子 (3)、抛掷一个不正规的骰子以剖析如何用蒙特卡罗方法模拟这些随机行为,以及基于MATLAB 软件的计算机实现。
5 抛掷一枚正规的硬币5.1 过程分析抛掷一枚硬币得到正面或反面的概率是1/2,我们可以把这种随机事件和[0,1]内的随机数建立联系。
概率是长期平均值,于是抛很多次时出现次数的比例接近0.5。
正面→[0,0.5],反面→(0.5,1]。
设x 为[0,1]内的随机数,f(x)定义如下:⎩⎨⎧≤<≤≤=15.05.00)(x x x f 反面,正面,f(x)将结果是正面或反面赋值到[0,1]内的一个数,随机赋值时我们可利用这个函随机数区间 出现的累积值 出现的百分比 x<0 0 0.00 0<x<0.5 0.5 0.50 5.2 符号说明n :抛掷硬币的次数COUNTER :记录得到正面的次数的计数器COUNTER/n :随机抛掷一枚正规的硬币得到正面向上的概率,n 的取值越大,得到的概率越接近理论值 i :第i 次抛掷硬币5.3 抛掷正规硬币的蒙特卡罗算法输入 模拟中生成的随机抛掷硬币的总次数n.. 输出 抛掷硬币时得到正面的概率. 第1步 初始化:COUNTER=0.第2步 对于i=1,2,…,n ,执行3,4步. 第3步 得到[0,1]内的随机数第4步 若5.00≤≤i x ,则COUNTER=COUNTER+1.否则,COUNTER 不变. 第5步 计算P (正面)=COUNTER/n. 第6步 输出正面的概率P (正面) 停止5.4MATLAB执行代码见附码(注意:重复执行代码,若不更改初始值,则每次执行,正面出现的次数会改变,正面出现的概率也随次数改变而改变)5.5结果分析下表给出了对于不同的n由随机数ix得到的结果:随着抛掷硬币次数n的增大,正面出现的概率也逐渐的接近0.5,即次数的一半。
当实验次数达到百万级以上时,模拟的概率值与理论值的误差仅为±0.0001。
6抛掷一个正规的骰子6.1过程分析一个骰子由点数{1,2,3,4,5,6}组成,抛掷一个正规的骰子,必须设计一种定义6个事件的方法。
而每个数值出现的可能性相等,所以每个数值出现的概率是1/6。
一个指定的数值出现的概率定义为:试验的总数数中指定的数值出现的次,6}{1,2,3,4,5将区间[0,1]等分为6个子区间,随机生成一个[0,1]上的随机数,它等可能的属于在6个子区间的某个区间上,这样就可将随机抛掷骰子出现的某个点数用随机生成[0,1]上的随机数属于某个子区间来模拟。
例如,随机生成的一个随机数属于[0,1/6],则认为随机抛掷一个正规的骰子出现的点数为1。
6.2符号说明n :抛掷骰子的次数COUNTER :记录次数的计数器COUNTER(j)/n :随机抛掷一个正规骰子的概率,n的取值越大,得到的概率越接近理论值i :第i次抛掷骰子j :骰子的点数6.3 抛掷一个正规骰子的蒙特卡罗算法输入 模拟中随机抛掷骰子的总数n. 输出 掷出{1,2,3,4,5,6}的百分比或概率.第1步 将COUNTER1到COUNTER6初始化为0. 第2步 对于i=1,2,…,n ,执行第3,4步 第3步 得到一个随机数,满足10≤≤i x .第4步 若i x 属于如下的区间,则相应的COUNTER+1.610≤≤i x COUNTER1=COUNTER1+1 6261≤<i x COUNTER2=COUNTER2+1 6362≤<i x COUNTER3=COUNTER3+1 6463≤<i x COUNTER4=COUNTER4+1 6564≤<i x COUNTER5=COUNTER5+1 165≤<i xCOUNTER6=COUNTER6+1第5步 计算掷出j={1,2,3,4,5,6}的概率为COUNTER(j)./n. 第6步 输出这些概率. 停止6.4 MATLAB 执行代码见附码(注意:重复执行代码,若不更改初始值,则每次执行,概率也会改变) 6.5 结果分析下表给出了掷10、100、1000、10000和100000次正规骰子的结果,掷1000005 0.3000 0.1700 0.1590 0.1758 0.1680 0.16676 0.1000 0.1400 0.1660 0.1682 0.1673 0.16677抛掷一个不正规的骰子7.1过程分析在抛掷一个正规的骰子的基础上考虑每个事件不是等可能出现的模型:抛掷一个不正规的骰子。
假定给骰子的某几个面加重使结果发生偏移。
将区间[0,1]分成6个不等长的子区间,子区间的长度代表骰子出现的某个值的概率。
例如,随机生成的一个随机数属于[0.4,0.7],对应着随机抛掷一个不正规骰子出现的点数为4的概率是0.3。
详细见下表:骰子出现的值出现的概率x的值赋值i1 0.1 [0.0,0.1] 12 0.1 (0.1,0.2] 23 0.2 (0.2,0.4] 34 0.3 (0.4,0.7] 45 0.2 (0.7,0.9] 57.2符号说明n :抛掷骰子的次数COUNTER :记录次数的计数器COUNTER(j)/n :随机抛掷一个不正规骰子的概率i :第i次抛掷骰子j :骰子的点数7.3抛掷一个不正规骰子的蒙特卡罗算法输入模拟中随机抛掷骰子的总数n.输出掷出{1,2,3,4,5,6}的百分比或概率.第1步将COUNTER1到COUNTER6初始化为0.第2步对于i=1,2,…,n,执行第3,4步第3步得到一个随机数,满足1≤x.0≤i第4步若x属于如下的区间,则相应的COUNTER+1.i≤x COUNTER1=COUNTER1+10≤1.0i<x COUNTER2=COUNTER2+11.0≤2.0ix COUNTER3=COUNTER3+1<2.0≤4.0i<x COUNTER4=COUNTER4+14.0≤7.0ix COUNTER5=COUNTER5+1<7.0≤9.0i<x COUNTER6=COUNTER6+19.0≤0.1i第5步计算掷出j={1,2,3,4,5,6}的概率为COUNTER(j)./n.第6步输出这些概率.停止7.4MATLAB执行代码见附码(注意:重复执行代码,若不更改初始值,则每次执行,概率也会改变)7.5结果分析当n=1000000,即重复进行100万次蒙特卡罗模拟试验后,得到的点数1,2,3,4,5,6的概率依次为p =0.0995 0.1002 0.1998 0.3002 0.2005 0.09998小结蒙特卡罗方法与一般计算方法有很大区别,一般计算方法对于解决多维或因素复杂的问题非常困难,而蒙特卡罗方法对于解决这方面的问题却比较简单。
其特点如下:直接追踪粒子,物理思路清晰,易于理解。
采用随机抽样的方法,较真切的模拟粒子输运的过程,反映了统计涨落的规律。
不受系统多维、多因素等复杂性的限制,是解决复杂系统粒子输运问题的好方法。
MC程序结构清晰简单。
研究人员采用MC方法编写程序来解决粒子输运问题,比较容易得到自己想得到的任意中间结果,应用灵活性强。
MC方法主要弱点是收敛速度较慢和误差的概率性质,其概率误差正比于,如果单纯以增大抽样粒子个数N来减小误差,就要增加很大的计算量。
蒙特卡罗算法已成为计算数学中不可缺少的组成部分,这主要是因为以下原因: 传统的分析方法受到了问题复杂性的限制。
MC方法直观,对实验者很有吸引力。
计算机变得更快更便宜。
量子理论的发展为我们提供了辐射与物质相互作用的截面数据。
9参考文献[1]叶其孝、姜启源,Frank R..Gioradano(美)等著,数学建模(原书第3版)[M],北京机械工业出版社,2005[2]MBA智库百科·蒙特卡罗方法/wiki/%E8%92%99%E7%89%B9%E5%8D%A1%E7%BD%9 7%E6%A8%A1%E6%8B%9F[3]百度文库·蒙特卡洛模拟法/view/4d5d6307cc175527072208a1.html [4]百度百科·蒙特卡罗模拟/view/284709.htm[5]新浪博客·Dwass (1957)蒙特卡罗模拟/s/blog_4b700c4c0100i28q.html10附码10.1抛掷一枚正规的硬币COUNTER=0;n=100; %硬币抛掷次数,赋初始值,大小可调节x=rand(1,n); %产生了n个随机数for i=1:nif x(i)>=0&&x(i)<=0.5COUNTER=COUNTER+1;endendCOUNTER %正面出现的次数P=COUNTER./n %正面出现的概率10.2抛掷一个正规的骰子for j=1:6COUNTER(j)=0;endn=10; %骰子抛掷次数,赋初始值,大小可调节x=rand(1,n); %产生了n个随机数for i=1:nif x(i)>=0&&x(i)<=1/6 COUNTER(1)=COUNTER(1)+1; elseif x(i)>1/6&&x(i)<=2/6 COUNTER(2)=COUNTER(2)+1; elseif x(i)>2/6&&x(i)<=3/6 COUNTER(3)=COUNTER(3)+1; elseif x(i)>3/6&&x(i)<=4/6 COUNTER(4)=COUNTER(4)+1; elseif x(i)>4/6&&x(i)<=5/6 COUNTER(5)=COUNTER(5)+1; elseif x(i)>5/6&&x(i)<=1 COUNTER(6)=COUNTER(6)+1; endendfor j=1:6p(j)=COUNTER(j)./n;endp %概率10.3抛掷一个不正规的骰子for j=1:6COUNTER(j)=0;endn=10; %骰子抛掷次数,赋初始值,大小可调节x=rand(1,n); %产生了n个随机数for i=1:nif x(i)>=0&&x(i)<=0.1 COUNTER(1)=COUNTER(1)+1; elseif x(i)>0.1&&x(i)<=0.2 COUNTER(2)=COUNTER(2)+1; elseif x(i)>0.2&&x(i)<=0.4 COUNTER(3)=COUNTER(3)+1; elseif x(i)>0.4&&x(i)<=0.7 COUNTER(4)=COUNTER(4)+1;elseif x(i)>0.7&&x(i)<=0.9 COUNTER(5)=COUNTER(5)+1; elseif x(i)>0.9&&x(i)<=1.0 COUNTER(6)=COUNTER(6)+1; endendfor j=1:6p(j)=COUNTER(j)./n;endp %概率。