蒙特卡罗方法学习总结

图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;

end

end 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误差

由中心极限定理可知，近似值与真值的误差为N

Z E Z N αλ<-)(?。式中的αλ的值可以根据给出的置信水平，查阅标准正态分布表来确定。

1.2.3收敛性与误差的关系

在一般情况下，求具有有限r 阶原点矩()∞

???

? ??--r

N 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 。此时伪随机数序列的最大容量2

2)(-=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 i

i P x F )(，)(x F 的直接抽样方法为)(,11

1∑∑=-=<≤=I i i I i i I F P P X X ξ。

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

3.2.2连续型分布的直接抽样方法

对于连续型分布，如果分布函数)(x F 的反函数)(1

x 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)(2

x 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 y

x xy +==θθ?，那么问题最终转换为在上半个单位圆内均匀抽样),(y x 的问题。为获得上半个单位圆内的均匀点，采用挑选法，在上半个单位圆的外切矩形内均匀投点。舍弃圆外的点，余下的就是所需要的点。

同样，借助MATLAB 程序来理解整个具体过程：

clear all;clc;

a=5.0;l=4.0;N=100000;

s=0;

for n=1:N

x=a*rand();

x=2*rand()-1;

y=rand();

while(x^2+y^2>1)

x=2*rand()-1;

y=rand();

end

sin_theta=(2*x*y)/(x^2+y^2);

if(x<=l*sin_theta)

s=s+1;

end

end

p=s*1.0/N;

result=2*l/(a*p);

fprintf('真值pi=%f\n',pi);

fprintf('计算pi=%f\n',result);

四、学习感想

在真实核物理实验中无可避免地涉及到放射性，为了避免对人员和环境造成伤害，降低成本，通常以模拟核物理实验来替代真实核物理实验。蒙特卡罗方法在实验核物理中的重要性不言而喻。

初次接触蒙特卡罗方法时，有很多基本的思路在之前都未曾接触过，所以会觉得蒙特卡罗方法非常难学。当然，本人在刚接触一门新的学科时经常会有这样的想法，一开始上课很难听得懂，作业也不会做。但我会尝试去接受、理解各种抽象的知识点，配合一定的练习去学习掌握好基本的知识。通过模仿老师所提供的或从其他途径获得的程序，慢慢地对蒙特卡罗方法有了一定的理解，上面的知识点是本人自认为是已经理解并掌握的。蒙特卡罗方法总的来说，可以总结出两个关键词：“随机抽样”、“统计”。

但是，有一些知识点（譬如蒙特卡罗方法在粒子运输中的应用），理解尚未足够透彻，难以形成自己的思路、组织语言写成学习总结。各学科之间经常会有着很多微妙的联系，在后续的学习中，将会接触到更多的专业知识，触类旁通，对理解蒙特卡罗方法应该会进一步拓宽和加深。

最后，非常感谢老师的谆谆教导！谢谢！