蒙特卡罗实验报告

蒙特卡罗方法
实验一
实验报告
蒙特卡罗方法实验一实验报告
一、实验目的
1、了解蒙特卡罗方法方法的基本思想；
2、掌握蒙特卡罗方法计算面积、体积的方法；
3、掌握由已知分布的随机抽样方法。

二、实验原理
Monte Carlo 方法，又称统计模拟方法或计算机随机模拟方法，是一种基于“随机数”进行数值模拟的方法，一种采用统计抽样理论近似求解物理或数学问题的方法。

倘若待求量可以表述成某些特征量的期望值、某些事件出现的概率或两者的函数形式，那么可采用蒙特卡罗方法求解。

在求解某些特征量的期望值或某些事件出现的概率时，必须构建合符实际的数学模型。

例如采用蒙特卡罗方法计算某函数所围面积时，构建的数学模型是构造一已知面积的可均匀抽样区域，在该区域投点，由伯努利定理大数定理可知，进入待求区域投点的频率依概率1收敛于该事件出现的概率（面积之比）。

由已知分布的随机抽样方法指的是由已知分布的总体中抽取简单子样。

抽样方法有： 直接抽样方法：离散型分布随机抽样方法、连续型分布直接抽样方法；挑选抽样方法;复合抽样方法；随机抽样一般方法：加抽样方法、减抽样方法、乘抽样方法、乘加抽样方法、乘减抽样方法、对称抽样方法、替换抽样方法、多为分布抽样方法、积分抽样方法；随机抽样其他方法：偏倚抽样方法、近似分布抽样方法、近似-修正抽样方法。

三、实验内容
1、安装所需计算工具（MA TLAB 、fortran 、C++等）；
2、编写一伪随机数发生器；（如乘加同余a=1366,c=150889,M=714025、a=9301,c=49297,M=233280；乘同余a=16807,M=232
-1；或采用其它方法）
以下内容选取一个采用自编伪随机数发生器进行计算，其余采用工具软件中自带伪随机数发生器进行计算。

3、求解以下区域的面积、体积：
3.1、给定曲线y =2 – x 2 和曲线y 3 = x 2，曲线的交点为：P 1( – 1，1 )、P 2( 1，1 )。

曲线围成平面有限区域，用蒙特卡罗方法计算区域面积；
3.2、计算22
22
11z x y z x y
⎧≥+⎪
⎨≤+--⎪⎩所围体积
其中{(,,)|11,11,02}x y z x y z Ω=-≤≤-≤≤≤≤。

4、对以下已知分布进行随机抽样： 4.1、()()
[]2
3
321,0,12
f x x x x =+
-∈；
4.2、()()
()[]11,1,21E f x f x x E k E =
⋅∈+
其中()()()()()2
123
221111211411ln 212221E x f x E x x x x E k E E E E E ⎧+-⎛⎫=+-+⎪ ⎪⋅⎝⎭⎪
⎨+⎡⎤⎪=-⋅+++-
⎢⎥⎪+⎣⎦⎩。

四、实验程序及其相关情况
第2题
function [SJS]=suiji(ZHONG) % SJS 产生的随机数
% ZHONG 输入的随机数种子 a=1366;c=150889;M=714025; SJS=mod(ZHONG*a+c,M)/M; 第3.1题
clear; clc;
M=0;%记录投点在所围图形中的个数 N= input('请输入总投点个数:\n'); ksi=0.89656; %用输入的方式 tic;
for i=1:N
ksi=suiji(ksi); x=2*ksi-1; ksi=suiji(ksi); y=2*ksi;
if y<=2-x^2
if y^3>=x^2 M=M+1; end end end toc
MIANJI=M/N*4 clear M N i x y;
计算结果：N =50000时面积为2.1431，计算时间约0.688s 。

第3.2题 clear;
clc;
M=0;%记录投点在所围图形中的个数 N= input('请输入总投点个数:\n'); tic;
for i=1:N
x=2*rand()-1;
y=2*rand()-1;
z=2*rand();
t=x^2+y^2;
s=z^2;
if s>=t
if t<=-s+2*z
M=M+1;
end
end
end
toc
MIANJI=M/N*8
clear M N i x y;
计算结果：N=50000时面积为3.1350，计算时间约0.282s。

第4.1题
clear;
clc;
M = input('输入所需产生随机变量的个数:\n');
x = zeros(M,1);
tic;
for i=1:M
if(rand()<=0.5)
x(i) = max(rand(),rand())
x(i) = max(x(i),rand());
x(i) = max(x(i),rand());
else
x(i) = min(rand(),rand());
x(i) = min(x(i),rand());
end
end
toc
clear M;
第4.2题
clear;
clc;
E = input('输入入射光子的能量(单位keV):\n');
E=E/511;
%计算系数
KE=(1-2*(E+1)/(E^2))*log(2*E+1)+0.5+4/E-1/(2*(2*E+1)^2); p1=(E/(2*E+1))^2/KE;
p2=p1+(2/E+2*E/(2*E+1)^2)/KE;
p3=p1+p2+(1-2*(E+1)/(E^2))*log(2*E+1)/KE;
tic;
a=rand(); if a<=p1
Cshe = (2*E+1)/(2*E*max(rand(),rand())-1); elseif a<=p2
Cshe = (2*E+1)/(2*E*rand()+1);
elseif a<=p3
Cshe = (2*E+1)^rand(); %可采用倒数抽样方法修正 else
Cshe = 2*E*rand()+1; end toc
clear E;
Cshe = Cshe*0.511
五、实验结果与分析
1、3.1和3.2抽样结果误差随抽样次数的关系图，及相关原因；
表1 实验记录表
序号
1
2 31 4
5
6
7
试验次数 103 5×104 5×104 1.2×105 1.5×105 1.8×106 2.0×107 试验时间 计算结果 实验误差
蒙特卡罗方法的近似值与真实值得误差问题，概率论中心极限定理指出，如果随机变量序列
Z 1，Z 2，。

，Z n 独立同分布，且具有有限异于零的方差，则
()dt e
21x lim x
x
-2
t
-
2
⎰
=⎪⎪⎭
⎫
⎝
⎛<-∧
π
σZ E N N
P Z
其中σ是随机变量Z 的均方差
()()()()
dt t f -t 2
2
2
）（⎰=
-=Z E Z E Z E σ
F （x ）是Z 的分布密度函数。

当N 充分大时，有如下近似式
απ
σ-1dt e
21x )(x
x
-2
t
-2
=≈⎪⎪⎭
⎫ ⎝
⎛<-⎰
∧
N Z E Z P N
其中α称为置信度，α-1成为置信水平。

该式的含义是：在α-1置信水平下，近似值与真值的误差为
N
σx 。

根据问题的要求确定出置信水平后，查标准正态分布表，确定出x ，就
可以得到近似值和真值的误差N
Z E Z N σx )(<
-∧
这样我们很容易发现误差的大小在σx 确定的情况下，N 越大误差越小
2、给出4.1和4.2的抽样框图、试验累积频率与理论累积频率关系图，并给出抽样次
数（>106
）与抽样时间。