蒙特卡罗方法学习总结

图1-1蒙特卡罗方法学习总结核工程与核技术2014级3班张振华20144530317一、蒙特卡罗方法概述1.1蒙特卡罗方法的基本思想1.1.1基本思想蒙特卡罗方的基本思想就是，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。

1.1.2计算机模拟打靶游戏为了能更为深刻地理解蒙特卡罗方法的基本思想，我们学习了蒲丰氏问题和打靶游戏两大经典例子。

下面主要对打靶游戏进行剖析、计算机模拟（MATLAB 程序）。

设某射击运动员的弹着点分布如表1-1所示，首先用一维数轴刻画出已知该运动员的弹着点的分布如图1-1所示。

研究打靶游戏，我们不用考察子弹的运动轨迹，只需研究每次“扣动扳机”后的子弹弹着点。

每一环数对应唯一确定的概率，且注意到概率分布函数有单调不减和归一化的性质。

首先我们产生一个在(0,1)上均匀分布的随机数（模拟扣动扳机），然后将该随机数代表的点投到P 轴上（模拟子弹射向靶上的一个确定点），得到对应的环数（即子弹的弹着点），模拟打靶完成。

反复进行N 次试验，统计出试验结果的样本均值。

样本均值应当等于数学期望值，但允许存在一定的偏差，即理论计算值应该约等于模拟试验结果。

clear all;clc;N=100000;s=0;for n=1:N %step 4.重复N 次打靶游戏试验x=rand();%step 1.产生在(0,1)上均匀分布的随机数if(x<=0.1)%step 2.若随机数落在(0.0,0.1)上，则代表弹着点在7环g=7;s=s+g;%step 3.统计总环数elseif(x<=0.2)%step 2.若随机数落在(0.1,0.2)上，则代表弹着点在8环g=8;s=s+g;elseif(x<=0.5)%step 2.若随机数落在(0.2,0.5)上，则代表弹着点在9环g=9;s=s+g;else%step 2.若随机数落在(0.5,1.0)上，则代表弹着点在10环g=10;s=s+g;endend gn_th=7*0.1+8*0.1+9*0.3+10*0.5;%step 5.计算、输出理论值fprintf('理论值：%f\n',gn_th);gn=s/N;%step 6.计算、输出试验结果fprintf('试验结果：%f\n',gn);1.2蒙特卡罗方法的收敛性与误差1.2.1收敛性由大数定律可知，应用蒙特卡罗方法求近似解，当随机变量Z 的简单子样数N 趋向于无穷大（N 充分大）时，其均值依概率收敛于它的数学期望。

1.2.2误差由中心极限定理可知，近似值与真值的误差为NZ E Z N αλ<-)(ˆ。

式中的αλ的值可以根据给出的置信水平，查阅标准正态分布表来确定。

1.2.3收敛性与误差的关系在一般情况下，求具有有限r 阶原点矩()∞<r Z E )21(<≤r 时的收敛速度为⎪⎪⎭⎫ ⎝⎛--rN O 1。

不难看出，在具有有限方差（2=r ）的假设下，试验次数N 需要增加两个数量级，所以在实际的应用中应作出适当的平衡。

1.3蒙特卡罗方法的优缺点1.3.1优点（1）可以逼真地模拟出具有随机性质的事物的特点及物理实验过程；（2）受几何条件的限制较小；（3）收敛速度与问题的维数无关，在处理多维度问题时有更明显的优势；（4）具有同时计算多个方案与多个未知量的能力；（5）误差容易确定；（6）程序结构简单，容易借助计算机实现。

1.3.2缺点（1）收敛速度慢，故通常不用来处理低维度问题；（2）误差具有概率性；（3）在粒子输运问题中，计算结果与系统大小有关。

综上所述，在应用蒙特卡罗方法解决现实问题时，应充分考虑蒙特卡罗方法的优缺点，尽可能地发挥其优点，同时还要注意到其缺点有可能带来的问题。

1.4蒙特卡罗方法在核科学技术方面的主要应用蒙特卡罗方法在核科学技术方面主要应用于粒子输运问题。

二、随机数2.1随机数的定义、性质及产生2.1.1随机数的定义由单位均匀分布（最简单、最基本的分布）抽取的简单子样称为随机数序列，其中每一个体称为随机数，通常用ξ表示。

2.1.2随机数的性质由随机数的定义，可以知道随机数本身具有独立性、均匀性两大性质。

2.1.3随机数的产生可以使用随机数表和物理方法来产生随机数。

随机数表，由0到9十个数字组成，每个数字等概率出现且相互独立。

若需要一个n 位有效数字的随机数，只需要将随机数表中每n 个相邻的随机数字合并在一起，且在最高位前加上小数点即可。

但是，随机数表在计算机中占用的内存空间很大，而且难以满足使用蒙特卡罗方法时需要产生大量随机数的要求，故该方法不适于在计算机上使用。

物理方法，利用某些物理现象，在计算机上增加某些特殊设备，可以在计算机上直接产生随机数。

但是，此方法产生的随机数序列无法重复实现，给结果验算带来极大的困难，而且须另外在计算机上联接附加设备，故该方法亦不适于在计算机上使用。

在现实使用蒙特卡罗方法时，通常使用数学方法来产生伪随机数来替代随机数。

2.2伪随机数2.2.1伪随机数在计算机上用递推公式 ,2,1),,,,(11==-+++n T k n n n k n ξξξξ产生随机数数列是最常见的产生伪随机数的方法。

对于给定的初始值k ξξξ,,,21 ，确定 ,2,1,=+n k n ξ。

2.2.2伪随机数的缺陷1）产生伪随机数的递推公式和处置确定后，整个随机数序列随之被唯一确定，不满足随机数相互独立的要求。

但是，只要递推公式设计比较合理，随机数间的相互独立性可以近似满足随机数的要求。

2）随机数序列是由递推公式确定得，而在计算机上所能表示的]1,0[上的数又是有限的。

如果出现有)(,n n n n ''<''''，使i n i n +''+'=ξξ，那么，随机数列就会出现周期性循环现象，不能满足随机数的要求。

我们可以定义伪随机数的周期为n n T '-''=，伪随机数的最大容量为n ''=λ。

但是，只要在一个循环周期的伪随机数个数多于所用随机数的个数时，此问题不会出现。

2.3产生伪随机数的乘同余方法2.3.1产生伪随机数的乘同余方法乘同余方法的一般形式：对于任一初始值1x ；伪随机数序列下面的递推公式确定：⎪⎩⎪⎨⎧==⋅≡+++ ,2,1,)(mod 111i M x M x a x i i i i ξ通常，取s M 2=（s 为计算机中二进制数的最大可能有效位数）；125+=k a （k 为计算机上所能容纳的最大整数）；11=x 。

此时伪随机数序列的最大容量22)(-=s M λ。

2.3.2C 语言乘同余方法产生伪随机数同样，我们可以应用C 语言来编程练习，加深对乘同余方法的理解。

摘取乘同余方法产生随机数模拟掷骰子代码中产生伪随机数的函数double randnum()及与之直接相关的代码：_int64M =pow(2,32);//step 1.给M 赋初始值M=2^s=2^32；_int64a =pow(5,13);//step 2.给a 赋初始值a=5^(2k+1)=5^13；_int64xi =1;//step 3.给x1赋值x1=1；double randnum(){_int64xn;double cauchy;{xn =(a*xi)%M;xi =xn;cauchy =xi*1.0/M;}//step 3.乘同余方法的递推公式；return cauchy;}三、由已知分布的随机抽样3.1随机抽样及其特点由已知分布的随机抽样指的是由己知分布的总体中抽取简单子样。

随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。

为方便起见，用F X 表示由己知分布函数)(x F 中产生的简单子样的个体。

对于连续型分布，常用分布密度函数)(x f 表示总体的己知分布，用f X 表示由己知分布密度函数)(x f 产生的简单子样的个体。

另外，在抽样过程中用到的伪随机数均称随机数。

3.2直接抽样方法3.2.1离散型分布的直接抽样方法对于任意离散型分布∑<=x x ii P x F )(，)(x F 的直接抽样方法为)(,111∑∑=-=<≤=I i i I i i I F P P X X ξ。

对于任意离散型分布，直接抽样方法是非常理想的。

在前面的蒙特卡罗方法基本思想中所详述的计算机模拟打靶游戏就是一个离散型分布的直接抽样的典例，在此不再赘述。

3.2.2连续型分布的直接抽样方法对于连续型分布，如果分布函数)(x F 的反函数)(1x F -存在，则直接抽样方法为：)(1ξ-=F X F 3.3替换法抽样3.3.1替换法抽样的定义为了实现某个复杂的随机变量y 的抽样，将其表示成若干个简单的随机变量n x x x ,,,21 的函数，),,,(21n x x x g y =得到n x x x ,,,21 的抽样后，即可确定y 的抽样，这种方法叫作替换法抽样。

3.3.2替换法抽样的应用能生动表示蒙特卡罗方法的基本思想的另一个典型问题——蒲丰氏投针试验求圆周率。

蒲丰氏投针试验在求解过程中需要用到θsin （θ在],0[π上均匀分布）的值，但在应用直接抽样时需要用到待求量π，我们考虑使用替换法抽样来解决这个矛盾。

ϕsin 服从分布：⎪⎩⎪⎨⎧≤≤--=其他,011,111)(2x x x f π直接抽样方法为：πξϕ2sin sin =，令θϕ2=，则θ在],0[π上均匀分布。

作变换⎩⎨⎧==θρθρsin cos y x （]1,0[∈ρ），则22sin y x y +=θ。

),(y x 表示上半单位圆内的点，若),(y x 在上半个单位圆内均匀分布，则对应地，θ在],0[π上均匀分布。

又222cos sin 2sin yx xy +==θθϕ，那么问题最终转换为在上半个单位圆内均匀抽样),(y x 的问题。

为获得上半个单位圆内的均匀点，采用挑选法，在上半个单位圆的外切矩形内均匀投点。

舍弃圆外的点，余下的就是所需要的点。

同样，借助MATLAB 程序来理解整个具体过程：clear all;clc;a=5.0;l=4.0;N=100000;s=0;for n=1:Nx=a*rand();x=2*rand()-1;y=rand();while(x^2+y^2>1)x=2*rand()-1;y=rand();endsin_theta=(2*x*y)/(x^2+y^2);if(x<=l*sin_theta)s=s+1;endendp=s*1.0/N;result=2*l/(a*p);fprintf('真值pi=%f\n',pi);fprintf('计算pi=%f\n',result);四、学习感想在真实核物理实验中无可避免地涉及到放射性，为了避免对人员和环境造成伤害，降低成本，通常以模拟核物理实验来替代真实核物理实验。