用蒙特卡洛方法估计积分方法及matlab编程实现

之杨若古兰创作
用蒙特卡洛方法估计积分
方法及matlab编程实现
专业班级：材料43
先生姓名：王宏辉
学号： 2140201060
指点教师：李耀武
完成时间： 2016年6月8日
用蒙特卡洛方法估计积分
方法及matlab编程实现实验内容：
1用蒙特卡洛方法估计积分
. 2用蒙特卡洛方法估计积分
.
请求：
（1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；
（2）利用计算机发生所选分布的随机数以估计积分值；
（3）进行反复试验，通过计算样本均值以评价估计的无偏性
；通过计算均方误差（针对第1类题）或样本方差（针对第2类题）
以评价估计结果的精度.
目的：
（1）能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等；
（2）熟练使用 MATLAB 对样本进行基本统计，从而获取数据的基本信息；
（3）能用 MATLAB 熟练进行样本的一元回归分析.
实验道理：
蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变量的数学期望，借助响应的随机数，利用样本均值估计数学期望，从而估计响应的积分值.
具体操纵如下：
X
;令Y=h(X)，则积分S=E（Y）；利用matlab软件，编程发生随机变量X的随机数，在由
Y的随机数，求出样本均值，以此估计积分值.
.
概率密度函数的拔取：
估计评价：
进行反复试验，通过计算样本均值以评价估计的无偏性；通过计
算均方误（针对第1类题，积得出）或样本方差（针对第2类题，积不出）以评价估计结果的精度.
程序设计：
根据成绩分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，响应程序设计成四类.
为了使程序具有普通性和方便当前使用：一重积分，程序保管为一个.m文本，被积函数，积分区间均采取键盘输入；二重积分，程序主体保管为一个.m文本，被积函数键盘输入，示性函数用function 语句构造，求分歧区域二重积分，只需改变function 函数内容.
编程完好解决用蒙特卡洛方法估计一重、二重积分值成绩.
程序代码及运转结果：
第一类一重积分程序代码：
%%%构造示性函数
function I=I1(x,a,b)
if x>=a&&x<=b
I=1;
else
I=0;
end
%%%%%%%%%%%%%%%%
%%第一类一重积分,程序主体：
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);%示性及求和向量
end
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]; %存放样本数字特征
运转结果：
%估计积分
1
m=f11
输入一元被积函数如x.*sin(x):x.*sin(x) g1 =
x.*sin(x)
输入积分下界a:0
输入积分上界b:pi/2
积分真值:1
输入样本容量 10^V1--10^V2:
V1:1
V2:5
n =
10
n =
100
n =
1000
n =
10000
n =
100000
Real =
1
EY1 =
D1 =
RD1 =
MSE1 =
m=
%估计积分
M=f11
输入一元被积函数如x.*sin(x):exp(-x.^2)
g1 =
exp(-x.^2)
输入积分下界a:0
输入积分上界b:+inf
输入样本容量 10^V1--10^V2: V1:1
V2:4
n =
10
n =
100
n =
1000
n =
10000
Real =
EY1 =
D1 =
RD1 =
MSE1 =
M =
第一类二重积分程序代码：
%%%构造示性函数，求分歧区域上积分只需更改示性函数
function I=I2(x,y)
if x^2+y^2<=1
I=1;
else
I=0;
end
%第一类二重积分程序主体
function outf12=f12()
g2=input('输入二元被积函数如exp(x.^2+y.^2):','s')%输入被积函数
g2=inline(g2,'x','y');
Real=input('积分真值:');%输入积分真值
fprintf('输入样本容量 10^V1*10^V1--10^V2*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);
y=randn(1,n);
for i=1:n
t2(i)=I2(x(i),y(i));%示性及求和向量
end
y2=g2(x,y)*(2*pi).*exp((x.^2+y.^2)/2);
Y2(j)=y2*t2'/n; %单次实验样本均值
end
t=ones(1,10);
EY=Y2*t'/10; %十次均值
D=abs(EY-Real); %绝对误差
RD=D/Real; %绝对误差
d=0;
for i=1:10
d=d+(Y2(i)-Real)^2;
end
d=d/(10-1);
EY2(m-V(1)+1)=EY; %样本容量为10^m时的样本均值
D2(m-V(1)+1)=D; %绝对误差
RD2(m-V(1)+1)=RD; %绝对误差
MSE2(m-V(1)+1)=d; %方差
end
Real,EY2,D2,RD2,MSE2
outf12=[EY2;D2;RD2;MSE2]; %存放样本数字特征
运转结果：
%估计积分
m=f12
输入二元被积函数如exp(x.^2+y.^2):exp(x.^2+y.^2)
g2 =
exp(x.^2+y.^2)
输入样本容量 10^V1*10^V1--10^V2*10^V2:
V1:1
V2:4
n =
10
n =
100
n =
1000 n =
10000 Real =
EY2 =
D2 =
RD2 =
MSE2 =
m =
第二类一重积分程序代码：
%%%构造示性函数
function I=I1(x,a,b)
if x>=a&&x<=b
I=1;
else
I=0;
end
%第二类一重积分程序主体
%
function outf21=f21()
g1=input('输入一元被积函数如exp(x.^2):','s')%输入被积函数
g1=inline(g1);
a=input('输入积分下界a:');%输入积分上上限
b=input('输入积分上界b:');
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);%示性及求和向量
end
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=0;
for i=1:10
d=d+(Y1(i)-EY)^2;
end
d=d/(10-1);
EY1(m-V(1)+1)=EY; %样本容量为10^m时的样本均值MSE1(m-V(1)+1)=d; %方差
end
EY1,MSE1
outf21=[EY1;MSE1]; %存放样本数字特征
%%%%
运转结果：
%估计积分
m=f21
输入一元被积函数如exp(x.^2):exp(x.^2) g1 =
exp(x.^2)
输入积分下界a:0
输入积分上界b:1
输入样本容量 10^V1--10^V2:
V1:1
V2:4
n =
10
n =
100
n =
1000
n =
10000
EY1 =
MSE1 =
m =
%用matlab 指令求积分
f=inline('exp(x.^2)')
f =
Inline function:
f(x) = exp(x.^2)
>> S=quadl(f,0,1)
S =
第二类二重积分程序代码：
%%%构造示性函数，求分歧区域上积分只需更改示性函数function I=I2(x,y)
if x^2+y^2<=1
I=1;
else
I=0;
end
%第二类二重积分函数主体
function outf22=f22()
g2=input('输入二元被积函数如 1./(1+x.^4+y.^4).^0.5:','s')%输入被积函数
g2=inline(g2,'x','y');
fprintf('输入样本容量 10^V1*10^V1--10^V2*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);
y=randn(1,n);
for i=1:n
t2(i)=I2(x(i),y(i));%示性及求和向量
end
y2=g2(x,y)*(2*pi).*exp((x.^2+y.^2)/2);
Y2(j)=y2*t2'/n; %单次实验样本均值
end
t=ones(1,10);
EY=Y2*t'/10; %十次均值
for i=1:10
d=d+(Y2(i)-EY)^2;
end
d=d/(10-1);
EY2(m-V(1)+1)=EY; %样本容量为10^m时的样本均值MSE2(m-V(1)+1)=d; %方差
end
EY2,MSE2
outf22=[EY2;MSE2]; %存放样本数字特征
运转结果：
%估计积分
m=f22
g2 =
输入样本容量 10^V1*10^V1--10^V2*10^V2:
V1:1
V2:4
n =
n =
100
n =
1000
n =
10000
EY2 =
MSE2 =
m =
实验结果清算：
第一类一重积分：
估计积分
积分真值：1积分估计值：
样本容量：10 100 1000 10000 100000
样本均值：
绝对误差：
绝对误差：
均方误差：
积分估计
样本容量：10 100 1000 10000
第一类二重积分：
积分估计值：
样本容量：10 100 1000 10000 样本均值：
绝对误差：
绝对误差：
均方误差：
第二类一重积分：
估计积分
积分估计 样本容量：10 100 1000 10000
第二类二重积分：
积分估计
样本容量：10 100 1000 10000
实验结果分析：
从第一类积分看，以估计积分
积分真值：1积分估计值： 样本容量：10 100 1000 10000 100000
样本均值：
绝对误差：
绝对误差：
均方误差：
随着样本容量的增大，样本均值有接近积分真值的趋势，绝对误差、绝对误差、均方误差呈减小趋势；随着样本容量的增大，样本均值有接近积分真值的趋势，说明估计具有没有偏性；绝对误差、绝对误差、均方误差呈减小趋势，说明增大样本容量能提高估计精度；验证了蒙特卡洛方法估计积分值的可行性，为后续估计第二类积分提供了参考.
从第二类积分看，以估计积分
样本容量：10 100 1000 10000
因为积分真值未知，没法直接比较估计值与积分值值；但随样本容量增大，样本方差减小，间接反映了估计精度的提高.蒙特卡洛方法估计值1.4590比拟用matlab
指令求得的积分结果1.4627，绝对偏差0.0038，绝对偏差0.0025.蒙特卡洛方法估计值与用matlab 指令求得的积分结果彼此验证.
总结与讨论：
蒙特卡洛方法是基于随机数的一种统计方法.
蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变量的数学期望，借助响应的随机数，利用样本均值估计数学期望，从而估计响应的积分值.
为使方法具有普通性，概率密度函数一重积分选择了
程序设计方面，本着使程序具有普通性和方便当前使用的准绳，根据成绩分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，响应程序设计成四类，并存储为.m文件，用蒙特卡洛方法估计积分值，一重积分只需调用响应程序即可；二重积分只需根据积分域点窜响应示性函数即可调用响应函数求解.
极大方便了同类成绩求解.
实验运转结果标明本方案可操纵性良好.
遗留成绩：本次实验未设计选用分歧概率密度函数，估计精度的比较，留有分歧条件下选用何种概率密度函数估计后果最好？如何缩短程序运转时间？如何对程序进行封装？如何更好评价第二类积分估计值无偏性和精度？等成绩.
姓名：王宏辉
班级：材料43
学号：2140201060。