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