用蒙特卡洛方法估计积分方法及matlab编程实现
- 格式:doc
- 大小:229.00 KB
- 文档页数:33
用蒙特卡洛方法估计积分方法及matlab编程实现
专业班级:材料43
学生姓名:王宏辉
学号:2140201060
指导教师:李耀武
完成时间:2016年6月8日
用蒙特卡洛方法估计积分 方法及matlab 编程实现
实验内容:
1用蒙特卡洛方法估计积分 2
0sin x xdx π⎰,2
-0x e dx +∞
⎰和
2
2
221
x
y x y e dxdy ++≤⎰⎰
的值,
并将估计值与真值进行比较。 2用蒙特卡洛方法估计积分 21
0x e dx ⎰
和
22x y +≤⎰⎰
的值,
并对误差进行估计。
要求:
(1)针对要估计的积分选择适当的概率分布设计蒙特卡洛方法;
(2)利用计算机产生所选分布的随机数以估计积分值; (3)进行重复试验,通过计算样本均值以评价估计的无偏性;通过计算均方误差(针
对第1类题)或样本方差(针对第2类
题)以评价估计结果的精度。 目的:
(1)能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数
及其期望、方差、协方差等;
(2) 熟练使用 MATLAB 对样本进行基本统计,从而获取数据的基本信息;
(3) 能用 MATLAB 熟练进行样本的一元回归分析。
实验原理:
蒙特卡洛方法估计积分值,总的思想是将积分改写为某个随机变量的数学期望,借助相应的随机数,利用样本均值估计数学期望,从而估计相应的积分值。 具体操作如下:
一般地,积分⎰=b
dx x g a )(S 改写成⎰⎰==b
b dx f h dx f g a
a )(x )(x )(x f(x))
(x S 的形
式,(其中为)f(x 一随机变量X 的概率密度函数,且)f(x 的支持域
)(}{b f ,a 0)(x |x ⊇>),f(x )
)
(x )(x g h =
);令Y=h(X),则积分S=E (Y );利用matlab 软件,编程产生随机变量X 的随机数,在由
⎩⎨
⎧∉∈==)
b (a,,)
b (a,,01I(x) ,)(x )(x y x x I h ,得到随机变量Y 的随机数,求出样本均值,以此估计积分值。
积分⎰⎰=A
dxdy g S )y (x,的求法与上述方法类似,在此不赘述。
概率密度函数的选取:
一重积分,由于要求)f(x 的支持域)(}{b f ,a 0)(x |x ⊇>,为使方法普
遍适用,考虑到标准正态分布概率密度函数2
2
e 21)(x x
f -=π
支持域为
R ,故选用2
2e 21)(x x f -
=
π
。
类似的,二重积分选用2
2
221)y (x,y x e f +-=π
,支持域为2R 。
估计评价:
进行重复试验,通过计算样本均值以评价估计的无偏性;通过计算均方误(针对第1类题,积得出)或样本方差(针对第2类题,积不出)以评价估计结果的精度。
程序设计:
依据问题分四类:第一类一重积分;第一类二重积分;第二类一重积分,第二类二重积分,相应程序设计成四类。
为了使程序具有一般性以及方便以后使用:一重积分,程序保存为一个.m文本,被积函数,积分区间均采用键盘输入;二重积分,程序主体保存为一个.m文本,被积函数键盘输入,示性函数用function 语句构造,求不同区域二重积分,只需改变function 函数内容。
编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。
程序代码及运行结果:
第一类一重积分程序代码:
%%%构造示性函数
function I=I1(x,a,b)
if x>=a&&x<=b
I=1;
else
I=0;
%保存为I1.m
%%%%%%%%%%%%%%%%
%%第一类一重积分,程序主体:
%保存为f11.m
function outf11=f11()
g1=input('输入一元被积函数如x.*sin(x):','s')%输入被积函数
g1=inline(g1);
a=input('输入积分下界a:');%输入积分上下限
b=input('输入积分上界b:');
Real=input('积分真值:');%输入积分真值
fprintf('输入样本容量10^V1--10^V2:\r')
V=zeros(1,2);
V(1)=input('V1:');%输入样本容量
V(2)=input('V2:');
for m=V(1):V(2)%样本容量10^m1--10^m2
n=10^m
for j=1:10
x=randn(1,n);
for i=1:n
t1(i)=I1(x(i),a,b);%示性及求和向量
y1=g1(x)*((pi*2)^0.5).*exp(x.^2/2);
Y1(j)=y1*t1'/n; %单次实验样本均值
end
t=ones(1,10);
EY=Y1*t'/10; %十次均值
D=abs(EY-Real); %绝对误差
RD=D/Real; %绝对误差
d=0;
for i=1:10
d=d+(Y1(i)-Real)^2;
end
d=d/(10-1);
EY1(m-V(1)+1)=EY; %样本容量为10^m时的样本均值
D1(m-V(1)+1)=D; %绝对误差
RD1(m-V(1)+1)=RD; %绝对误差
MSE1(m-V(1)+1)=d; %方差
end
Real,EY1,D1,RD1,MSE1
outf11=[EY1;D1;RD1;MSE1]; %存放样本数字特征
%保存为f11.m