《数理统计》课程设计题目：用M-C 方法求积分1() f x dx⎰【题目要求：f（x)自定，n≥500,考虑n对结果的影响，即做多组n下的模拟值，并作模拟值与n的散点图，同时比较模拟值与真实值的差异，散点图表示。
并做差异值序列的描述性统计（均值、方差、标准差、峰度系数、偏度系数、众数、中位数、四分位数等）。
积分区间可根据需要调整。
】学院：数学学院专业班级：应用数学09-2班姓名：李明学号： 20096312指导教师：谭常春2012.6.20一、M-C方法概述M-C方法即蒙特卡洛方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。
这一方法源于美国在第二次世界大战中研制原子弹的“曼哈顿计划”。
该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。
该方法基本思想很早以前就被人们所发现和利用。
17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。
19世纪人们用投针试验的方法来决定π。
高速计算机的出现，使得用数学方法在计算机上大量模拟这样的试验成为可能。
其实质是通过大量随机试验，利用概率论解决问题的一种数值方法，基本思想是基于概率和体积间的相似。
Monte Carlo方法计算结果收敛的理论依据来自于大数定律，且结果渐进地服从正态分布的理论依据是中心极限定理。
以上两个属性都是渐进性质，要进行很多次抽样，此属性才会比较好地显示出来，如果Monte Carlo计算结果的某些高阶距存在，即使抽样数量不太多，这些渐进属性也可以很快地达到。
二、M-C方法与数值积分用数值积分方法计算积分，如21()xx f x dx⎰，如果我们能够得到f(x)的原函数F(x)，那么直接由表达式: F(x2)-F(x1)可以得到该定积分的值。
但是，很多情况下，由于f(x)太复杂，无法计算得到原函数F(x)的显式解，这时我们就只能用数值积分的办法。
数值积分的基本原理是在自变量x的区间上取多个离散的点，用单个点的值来代替该小段上函数f(x)值。
常规的数值积分方法是在分段之后，将所有的矩形小块的面积全部加起来，用这个面积来近似函数f(x)与x轴围成的面积。
这样做当然是不精确的，但是随着分段数量增加，误差将减小，近似面积将逐渐逼近真实的面积。
Monte Carlo方法和上述类似。
差别在于，Monte Carlo方法中，我们不需要将所有方柱的面积相加，而只需要随机地抽取一些函数值，将他们的面积累加后计算平均值就够了。
随着抽取点增加，近似面积也将逼近真实面积。
三、M-C方法的形式与一般步骤做Monte Carlo时，求解积分的一般形式是：21()()xx f x x dψ⎰；x为自变量，它应该是随机的，定义域为(x1, x2)，f(x)为被积函数，ψ(x)是x的概率密度。
Monte Carlo方法分为一下四个个步骤：1．依据概率分布ψ(x)不断生成随机数x, 并计算f(x)：由于随机数性质，每次生成的x 的值都是不确定的，为区分起见，我们可以给生成的x 赋予下标。
如x i 表示生成的第i 个x 。
生成了多少个x ，就可以计算出多少个f(x)的值。
2．将这些f(x)的值累加，并求平均值例如我们共生成了N 个x ，这个步骤用数学式子表达就是3．到达停止条件后退出。
常用的停止条件有两种：一种是设定最多生成N 个x ，数量 达到后即退出；另一种是检测计算结果与真实结果之间的误差，当这一误差小到某个范围之内时退出。
4．误差分析：Monte Carlo 方法得到的结果是随机变量，因此，在给出点估计后，还需要给出此估计值的波动程度及区间估计。
严格的误差分析。
首先要从证明收敛性出发，再计算理论方差，最后用样本方差来替代理论方差。
四、常见随机数的生成及相关函数1、rand() 生成（0,1）区间上均匀分布的随机变量。
基本语法：rand([M,N,P ...]) 生成排列成M*N*P... 多维向量的随机数。
如果只写M ，则生成M*M 矩阵；如果参数为[M,N]可以省略掉方括号。
2、randn() 生成服从标准正态分布（均值为0，方差为1）的随机数。
基本语法和rand()类似：randn([M,N,P ...]);生成排列成M*N*P... 多维向量的随机数。
如果只写M ，则生成M*M 矩阵；如果参数为[M,N]可以省略掉方括号。
3、unifrnd() 和rand()类似，rand （）可以看作其特殊情况。
这个函数生成某个区间内均匀分布的随机数。
基本语法：unifrnd(a,b,[M,N,P,...]);生成的随机数区间在(a,b)内，排列成M*N*P... 多维向量。
如果只写M ，则生成M*M 矩阵；如果参数为[M,N]可以省略掉方括号。
4、normrnd() 和randn()类似，randn （）可以看作其特殊情况。
此函数生成指定均值、标准差的正态分布的随机数。
基本语法：normrnd(mu,sigma,[M,N,P,...])；生成的随机数服从均值为mu ，标准差为sigma （注意标准差是正数）正态分布，这些随机数排列成M*N*P... 多维向量。
如果只写M ，则生成M*M 矩阵；如果参数为[M,N]可以省略掉方括号。
5、chi2rnd() 此函数生成服从卡方分布的随机数。
卡方分布只有一个参数：自由度v 。
基本语法：chi2rnd(v,[M,N,P,...])；生成的随机数服从自由度为v 的卡方分布，这些随机数排列成M*N*P... 多维向量。
如果只写M ，则生成M*M 矩阵；如果参数为[M,N]可以省略掉方括号。
6、frnd() 此函数生成服从F 分布的随机数。
F 分布有2个参数：v1, v2。
基本语法：frnd(v1,v2,[M,N,P,...])；生成的随机数服从参数为(v1,v2)的卡方分布，这些随机数排列成M*N*P... 多维向量。
如果只写M ，则生成M*M 矩阵；如果参数为[M,N]可以省略掉方括号。
1()Nii f x N=∑7、trnd() 此函数生成服从t 分布，t 分布有1个参数：自由度v 。
基本语法：trnd(v,[M,N,P,...])；生成的随机数服从参数为v 的t 分布，这些随机数排列成M*N*P... 多维向量。
如果只写M ，则生成M*M 矩阵；如果参数为[M,N]可以省略掉方括号。
8、其他常见分布：生成beta 分布随机数的语法是：betarnd(A,B,[M,N,P,...]) 生成Gamma 分布随机数的语法是：gamrnd(A,B,[M,N,P,...]) 生成指数分布随机数的语法是：betarnd(mu,[M,N,P,...])。
五、积分的计算计算定积分： 20x e dx⎰（1）数学方法：我们已知xe 的原函数是xe ，那么定积分值就是：2e e -=6.3890561 。
计算这个数值可以在Matlab 中输入代码：exp(2)-exp(0)。
得到的值是此定积分的真实值。
（2）常规数值积分：在(0,2)x ∈区间内取N 个点，计算各个点上的函数值，然后用函数值乘以每个区间宽度，最后相加。
Matlab 代码：N=100;x=linspace(0,2,N); sum(exp(x).*2/N)试着调大N 的值，最后的结果将更接近真实值。
（3）Monte Carlo 积分法：在(0,2)x ∈内随机取N 个点，计算各个点上的函数值，最后求这些函数值的平均值再乘以2。
看Matlab 代码：N=100;x=unifrnd(0,2,N,1); mean(2*exp(x))同样的，通过增大N ，这种方法得到的结果也将越来越接近真实值。
对M-C 方法的解释：这个积分要求的积分形式是：20x e dx⎰， 还不完全是21()()x x f x x dx ψ⎰形式，故先做变换，21(2)()2xe dx ⎰，这里2xe 是f(x)；1/2是ψ(x)，它表示，在取值范围(0,2)区间内，x 服从均匀分布。
显然若是积分区间为（0，1),则是在区间（0，1）内服从均匀分布，即从0到1内随机的取出N 个点。
对于其代码语句解释如下：N=100 设定停止条件，共做N 次Monte Carlo 模拟；x=unifrnd(0,2,N,1) 按照(0,2)区间均匀分布概率密度对x随机抽样，共抽取N 个x。
mean(2*exp(x)) 2*exp(x)作用是对每个x i计算f(x i)的值,共可得到N个值，这个相当于第一个步骤后半部分；Mean()函数的作用是将所有的f(x i)加起来取平均值。
六、M-C方法对积分计算的分析对上述函数计算中，N的值从500至1000每隔10取定一个值（N值即区间0到2内取的x点数），分别计算其积分模拟值。
并作模拟值与N的散点图，同时比较模拟值与真实值的差异，散点图表示。
并做差异值序列的描述性统计（均值、方差、标准差、峰度系数、偏度系数、众数、中位数、四分位数等）。
通过计算得出下列结果：第一列为N的值，第二列为对应的模拟值，第三列为N对应的模拟值与真实值的差异。
（真实值为6.3891）500 6.0793 0.30976 510 6.2867 0.10232 520 6.6887 -0.29968 530 6.2313 0.15771 540 6.4241 -0.03501 550 6.3359 0.053158 560 6.5737 -0.18465 570 6.5098 -0.12075 580 6.2231 0.16598 590 6.5384 -0.14935 600 6.3856 0.0034121 610 6.4427 -0.053672 620 6.4766 -0.087552 630 6.2921 0.096933 640 6.2444 0.14465 650 6.4517 -0.06269 660 6.4666 -0.077554 670 6.3417 0.047389 680 6.4812 -0.092182 690 6.5255 -0.13645 700 6.3259 0.063167 710 6.5512 -0.16213 720 6.4639 -0.074884 730 6.4253 -0.036203 740 6.5366 -0.1475 750 6.5665 -0.17745 760 6.5481 -0.15909 770 6.5607 -0.17165 780 6.3552 0.033831 790 6.4636 -0.074584 800 6.3712 0.017859 810 6.4532 -0.064155 820 6.5182 -0.12912 830 6.1661 0.22291 840 6.1409 0.24814 850 6.4171 -0.028075 860 6.2438 0.1453 870 6.4517 -0.062662 880 6.266 0.12309 890 6.2639 0.12511 900 6.4105 -0.021448 910 6.2054 0.18363 920 6.5427 -0.15366 930 6.3731 0.01594 940 6.5718 -0.18274 950 6.2956 0.093458 960 6.3882 0.00088631 970 6.3696 0.019428 980 6.4955 -0.10648 990 6.3351 0.053931 1000 6.3671 0.021975计算第三列差异值的描述性统计可得：均值：-0.032097 方差：0.018636 标准差：0.13651 中位数: -0.042936模拟值与N（N值即区间0到2内取的x点数）的散点图：（真实值为6.3891）差异值与N的散点图显然，此时由于点数取的较少，收敛性并不明显，下面是取点较多的情况：从500至5500每隔10取一个N值，共501个N值：从500至10000每隔10取一个N 值，共951个N值：从500至20000每隔10取一个N 值，共1951个N值：从上面的图中可以发现，N 值越大，即积分区间内所取的点数越多，模拟值越接近真实值，且模拟值的波动区间也越小。