概率论与数理统计实验报告
——蒙特卡洛算法计算积分
实验时间:
姓名:司默涵
学号:2110505018
班级:计算机11
一、实验目的
(1)能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数及其
期望、方差、协方差等;
(2) 熟练使用 MATLAB 对样本进行基本统计,从而获取数据的基本信息;
(3) 能用 MATLAB 熟练进行样本的一元回归分析。
二、实验要求
(1)针对要估计的积分选择适当的概率分布设计蒙特卡洛方法;
(2)利用计算机产生所选分布的随机数以估计积分值;
(3)进行重复试验,通过计算样本均值以评价估计的无偏性;通过计算均方误
差(针对第1类题)或样本方差(针对第2类题)以评价估计结果的精度。
三、实验原理
1. 蒙特卡洛法的思想简述
当我们所求解问题是某种随机事件出现的概率,或者是某个随机变量的期望值时,通过某种“实验”的方法,以这种事件出现的频率估计这一随机事件的概率,或者得到这个随机变量的某些数字特征,并将其作为问题的解。 有一个例子我们可以比较直观地了解蒙特卡洛方法:假设我们要计算一个不规则图形的面积,那么图形的不规则程度和分析性计算(比如,积分)的复杂程度是成正比的。蒙特卡洛方法是如下计算的:假想有一袋豆子,把豆子均匀地朝这个图形上撒,然后数这个图形之中有多少颗豆子,这个豆子的数目就是图形的面积。当豆子越小,撒的越多的时候,结果就越精确。在这里我们要假定豆子都在一个平面上,相互之间没有重叠。
2. 蒙特卡洛法与积分
通常蒙特卡洛方法通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题,蒙特卡洛方法是一种有效的求出数值解的方法。一般蒙特卡洛方法在数学中最常见的应用就是蒙特卡洛积分。非权重蒙特卡洛积分,也称确定性抽样,是对被积函数变量区间进行随机均匀抽样,然后对被抽样点的函数值求平均,从而可以得到函数积分的近似值。此种方法的正确性是基于概率论的中心极限定理。
3. 本实验原理简述
在本实验中,我们主要是计算积分值与误差比较。在计算积分时,我们要选择合适的变量分布,其中有均匀分布,有正态分布,要视情况而选择。在利用蒙特卡洛方法计算积分时,我们要分情况。 ①对于积分为这种形式,我们可以转化为
这种形式,然后利用其等于(b-a)E(x)的计算结果。E(x)可利用求随机变量的均值来得到。
②对于积分为这种形式的,我们依然可以利用上面的算法计算,将其化为
y)可利用求随机变量的均值来得到。
?对于有无穷的区间,我们应选择正态分布。 即可得。E(x,
三.实验内容和实验结果(附图)
1.计算下面的积分,并与真值比较 ?
2??
①?xsinxdx ②
0?e?xdx2 ③.0x?y?12??ex2?y2dxdy2
2.计算下面的积分,并与平均值比较计算方差
1
④. ?e
0x2dx ⑤.
x2?y2???1dxdy
计算过程及程序如下 ?
2
①
分析:首先产生若干组按在[0,π/2]上的均匀分布的简单随机样本,然后计算样本对应的y=xsinx的值,求出y的均值,最后用样本的均值估计随机变量的期望。即得到一个积分值的估计,求出若干组估计值根据其计算出估计值的均值,得到最后的估计。用此估计和正确值比较。
n=10000;m=10;S=0;I=0;
r=unifrnd(0,pi/2,m,n)
for j=1:m
s=0;
for i=1:n
s=s+r(j,i)*sin(r(j,i));
end ?xsin0xdx
S(j)=pi/2*s/n
end
D=0;d=0;
for j=1:m
D=D+(S(j)-1)^2;
end
d=D/(m-1)
结果S =1.009 1.0023 0.9898 0.9859 0.9948 0.9828 0.9963 1.0049 0.9943 0.9924
d =8.5885e-005
??
②
分析:题目中因为积分区间无限故不能使用均匀分布,此时可以利用matlab产生正态分布模拟。此次产生1000组每组100个的样本,观察此次的误差及方差。
x=randn(100,1000); %产生满足正态分布的100组每组1000个样本的简单随机样本
x_aver=sum(x,2)/1000;
e=exp(-x_aver.^2);
ans=sqrt(pi/2)*e;
plot(ans)
ans_aver=sum(ans,1)/100
corans=sqrt(pi/2)
wucha=abs(ans_aver-corans)
fangcha=(ans-corans).^2; ?e?x2dx0
fangcha=sum(fangcha,1)/99
结果ans_aver =1.2520 corans = 1.2533 wucha = 0.0013 fangcha =0.0020 ③.x2?y2?1??ex2?y2dxdy
m=10000;sum=0;n=50;D=0;
X=unifrnd(-1,1,n,m);Y=unifrnd(-1,1,n,m);
for i=1:n
a=0;
for j=1:m
if(X(i,j)^2+Y(i,j)^2<=1)
Z(i,j)=exp(X(i,j)^2+Y(i,j)^2);
a=a+Z(i,j);
end
end
S(i)=a/m;sum=sum+S(i); end
I=sum/n*4
for i=1:n
D=D+(S(i)*4-pi*(exp(1)-1))^2; end
d=D/n
结果I =5.4089;d =0.0011
1
④ e?.
0x2dx
x=rand(100,1000);
y=exp(x.^2);
y_aver=sum(y,2)/1000;
ans=sum(y_aver,1)/100 fangcha=(y_aver-ans).^2; fangcha=sum(fangcha,1)/99 plot(y_aver)
结果ans =1.4627;fangcha = 0.0171
⑤.
x2?y2???1dxdy
n=1000;m=100;sum=0;S=0;I=0;
x=unifrnd(-2,2,m,n);y=unifrnd(-2,2,m,n); for j=1:m
s=0;
for i=1:n
if x(j,i)^2+y(j,i)^2<=4
s=s+16/sqrt(1+x(j,i)^4+y(j,i)^2); end
end
S(j)=s/n;sum=sum+S(j);
end
I=sum/m;D=0;d=0;
for j=1:m
D=D+(I-S(j))^2;
end
d=D/(m-1)
结果I =7.4531;d =0.0188
四.实验心得
(1)由实验结果可以看出,用蒙特卡罗估计算法进行积分是可行并且有效的,
可以用它来计算原函数难以求得的被积函数的积分及广义积分的值。与其他用来计算积分的方法如simpson、newton等相比它具有简单易懂且容易实现的特点。
b
(2)从理论上讲,给一个积分?g(x)*f(x)dx,我们可以用很多种分布去模拟计
a
算,只需要选取的分布拥有概率密度g(x)即可。如在本次实验中,计算?
?xe?dx时,可用正态分布模拟,此时g(x)?
02x?1?2e,f(x)?2?*e2,21x22
也可用均匀分布进行模拟,此时g(x)?1,f(x)?e?x21?x?2e。两种模拟x
计算的结果差不多,前者为0.888373,后者为0.86197,与真值2?0.88623很接近。
(3)对第一和第二两个积分的模拟结果求方差,可以看出两者都在10?6数量级,
这说明用蒙特卡罗模拟估算法的可靠性很高,与真值极为接近。