用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx

1、（3） 进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针 对第1类题）或样本方差（针对第2类题）以评价估计结果的精度。目的：（1）能通过MATLAB或其他数学软件了解随机变量的概率密度、分布函数 及其期望、方差、协方差等；（2）熟练使用MATLAB对样本进行基本统计，从而获取数据的基本信息；（3）能用MATLAB熟练进行样本的一元回归分析实验原理：蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变 量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从 而估计相应的积分值。具体操作如下:式，（其中为f（x） 一随机变量 X的概率密度函数，且f（x）的支持域x|

2、f（x）O=（ a,b） , h（x-gx）；令 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（xF彳 廿类似的，二重积分选用f（x, y） e 2，支持域为R22兀估计评价:进行重复试验，通过计算样本均值以评价估计的无偏性；通过计 算均方误（针对第

3、1类题，积得出）或样本方差（针对第2类题，积 不出）以评价估计结果的精度。程序设计：依据问题分四类：第一类一重积分；第一类二重积分；第二类一 重积分，第二类二重积分，相应程序设计成四类。为了使程序具有一般性以及方便以后使用：一重积分，程序保存 为一个.m文本，被积函数，积分区间均采用键盘输入；二重积分， 程序主体保存为一个.m文本，被积函数键盘输入，示性函数用 fun ction 语句构造，求不同区域二重积分，只需改变 fun ction 函数内容。编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。程序代码及运行结果：第一类一重积分程序代码：%构造示性函数fun ctio n I=l1（x,

4、a,b）if x=a&x=bl=1;elsel=0;end%保存为I1.m%第一类一重积分，程序主体：%保存为f11.mfun ctio n outf11=f11（）输入被积g1=input（输入一元被积函数如 x.*sin（x）:,s）% 函数g1=i nlin e（g1）;a=input（输入积分下界a:）;%输入积分上下限b=input（输入积分上界b:Real二input（积分真值:%输入积分真值 fprintf（ 输入样本容量 10AV1-10AV2:r） V=zeros（1,2）;V（1）=i nput（V1:% 输入样本容量V（2）=i nput（V2:for m=V（1）:V（

5、2）% 样本容量 10八口1-10八口2n=10八口for j=1:10x=ra ndn （1, n）;for i=1: nt1（i）=l1（x（i）,a,b）;% 示性及求和向量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;d=d+（Y 1（i）-Real）A2;d=d/（10-1）;EY 1（m-V（1）+1）=E Y; % 样本容量为10八m时的样本均值D1（m-V（1）

6、+1）=D;RD1（m-V（1）+1）=RD;MSE1（m-V（1）+1）=d; % 方差Real,E Y1,D1,RD1,MSE1outf11=EY1;D1;RD1;MSE1; % 存放样本数字特征运行结果：2%古计积分 xs in xdx，积分真值为1m=f11输入一元被积函数如x.*sin（x）:x.*sin（x）g1 = x.*si n（x）pi/2输入样本容量10AV1-10AV2:5n 二100100010000100000Real =EY1 =1.2635 1.0088 1.0066 1.0109 1.0018D1 =0.26350.00880.00660.01090.0018R

7、D1 =MSE1 =0.64390.02050.00280.00060.0001m=1.26351.00881.00661.01091.0018%估计积分 e-x2dx真值为0.8862M=f11输入一元被积函数如x.*si n（x）:exp（-x42）g1 =exp（-x.A2）+inf积分真值：piA0.5/2%0.886240.88620.93330.90770.88730.88710.04700.02150.00100.00090.05310.02430.00120.19270.01120.00160.0000M =0.1927 0.0112 0.0016第一类二重积分程序代码:%造示

8、性函数，求不同区域上积分只需更改示性函数fun cti on I=l2（x,y）if xA2+yA2=1%保存为I2.m%第一类二重积分程序主体%保存为f12.mfun ctio n outf12二f12（）g2=in put（输入二元被积函数如 exp（x.A2+y.A2）:）%入被积函数g2=i nlin e（g2,xy%输入积分真值fprintf（ 输入样本容量 10AV1*10AV1-10AV2*10AV2:）V=zeros（1,2）;y=ra ndn （1, n）;t2（i）=l2（x（i）,y（i）;y2=g2（x,y）*（2*pi）.*exp（x.A2+y.A2）/2）;Y 2（

9、j）=y2*t2t=o nes（1,10）;EY二Y 2*td=d+（ Y2（i）-RealF2;EY 2（m-V（1）+1）=E Y;%样本容量为10a m时的样本均值D2（m-V（1）+1）=D;绝对误差RD2（m-V（1）+1）=RD;MSE2（m-V（1）+1）=d;方差Real,E Y2,D2,RD2,MSE2outf12=EY2;D2;RD2;MSE2;%古计积分 ex2 y2dxdy，真值为 pi*（exp（1）-1）%5.3981xy2m=f12exp（x.A2+y.A2）g2 =pi*（exp（1）-1）%5.3981输入样本容量 10AV1*10AV1-10AV2*10AV2:5.3981EY2 =4.77025.12505.43175.4041D2 =0.62790.27320.03350.0060RD2 =0.11630.05060.00620.0011MSE2 =3.89650.55640.02470.0017m =第二类一重积分程序代码:I=1;I=0;%第二类一重积分程序主体%程序保存为f21.mfun