统计计算统计模拟π的估计

1、统计模拟作业班级学号姓名作业一：的估计1. 分析方法 从十七世纪中叶起，人们开始用更先进的分析方法来求的近似值，其中应用的主要工具是收敛的无穷乘积和无穷级数，下面列举一些用此类方法求近似值的实例。(1)由公式 推出=4编写程序：syms kx=symsum(-1)k/(2*k+1),k,0,10)y=4*x得出当k=10时，3.232315809405593编写程序syms kx=symsum(-1)k/(2*k+1),k,0,20)y=4*x得出当k=20时，3.189184782277595依次，加大k的值K=50，3.161198612987050K=100，3.151493401070

2、990K=200，3.146567747182986K103.232315809405593503.1611986129870501003.1514934010709902003.146567747182986（2）沃里斯(Wallis)方法在积分学中我们经常会遇到如下的沃利斯(Wallis)公式。沃利斯(Wallis)公式揭示了与整数之间的一种很不寻常的关系。但在实际学习中很少注意到沃利斯(Wallis)公式,更不会关注它的应用。实际上,沃利斯(Wallis)公式有许多作用,经常有以下几方面的应用。1应用于极限计算中。由于沃利斯(Wallis)公式与极限有关,所以有些极限的计算可以通过沃利斯

3、(Wallis)公式很容易计算出来。2.应用于积分计算中。对于一些用积分法不易求得原函数的积分，而用沃利斯(Wallis)公式却很容易解决问题。编写程序：format long x=1;for k=1:10 x=x*(2*k/(2*k-1)*2*k/(2*k+1);endy=2*x得k=10时，3.067703806643498增加k的值K=20，3.103516961539230K=50，3.126078900215409K=100，3.133787490628159K=10000，3.141514118681864K=1000000，3.141591868191880k103.067703

4、806643498203.103516961539230503.1260789002154091003.133787490628159100003.14151411868186410000003.141591868191880（3）蒙特卡罗法 取一正方形A，以A的一个顶点为圆心，A的边长为半径画圆，取四分之一圆（正方形内的四分之一圆）为扇形B。已知A的面积，只要求出B的面积与A 的面积之比，就能得出，再由B的面积为圆面积的四分之一，利用公式即可求出的值。因此，我们的目的就是要找出的值。可以把A和B看成是由无限多个点组成，而B内的所有点都在A内。随机产生个点，若落在B内的有个点（假定A的边长为1

5、，以扇形圆心为坐标系原点。则只要使随机产生横纵坐标、满足的点，就是落在B内的点），则可近似得出的值，即，由此就可以求出的值。编写程序：i=1;m=0;n=1000;for i=1:n a=rand(1,2); if a(1)2+a(2)2<=1 m=m+1; endendp=vpa(4*m/n,30)程序运行结果：p =3.14000000000000000000000000000（4） 泰勒级数法函数的泰勒级数展开式为：将代入上式有.利用这个式子就可以求出的值了。编写程序：i=1;n=1000;s=0;for i=1:n s=s+(-1)(i-1)/(2*i-1);endp=vpa(4

6、*s,30)程序运行结果：p =3.14059265383979413499559996126当取的值为10000时，就会花费很长时间，而且精度也不是很高。原因是时，的展开式收敛太慢。因此就需要找出一个使得收敛更快。若取，则我们只有找出与的关系，才能求出的值。令，根据公式有，则有。所以可以用来计算的值。（5）利用公式 推出=4()编写程序：i=1;n=1000;s=0;s1=0;s2=0;for i=1:n s1=s1+(-1)(i-1)*(1/2)(2*i-1)/(2*i-1); s2=s2+(-1)(i-1)*(1/3)(2*i-1)/(2*i-1);ends=s1+s2;p=vpa(4*

7、s,30)程序运行结果：p =3.14159265358979323846264338328 显然，级数收敛越快，取同样的值可以得到更高的精度。以同样的方法，能得出，程序和上面的一样。这样的近似值可以精确到几百位。（6） 麦琴公式由推出=4()，这个公式由英国天文学教授约翰·麦琴于1706年发现。他利用这个公式计算到了100位的圆周率。麦琴公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。编写程序：syms n;f1=(-1)(n-1)*(1/5)(2*n-1)/(2*n-1);f2=(-1)(n-1)*(

8、1/239)(2*n-1)/(2*n-1);ans1=symsum(f1,n,1,28);ans2=symsum(f2,n,1,28);ans=vpa(4*(4*ans1-ans2),100)得3.141592653589793238462643383279502884197130451046268578972203255663716036677133432949011735665451127（7）用下列公式 求值，并了解的其他公式。编写程序：digits(10000);N=100000000;p=1;for m=1:N;m1=m*m;m2=2*m;n=(4*m1)/(m2-1)*(m2+1)

9、;p=p*n;endp=vpa(p*2);改变N的值，分析得到的不同值的大小与N的关系。PI=3. 141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825 当N=10*104时，=3.14151411868195662435709891724400222301483154296875当N=10*108时，=3.141592643066262180440162410377524793148040771484375 在使用MATLAB的过程中发现，当N值越大的时候，PI的值

10、精确度越高，然而N大到一定范围后得出结果的过程也就变得非常慢，并且在10*108后每增加10倍后PI没有很大的变化，并且运算速度也变得非常慢了。2.概率方法 取一个二维数组 (x,y) ,取一个充分大的正整数n,重复n次，每次独立的从（0,1）中随机地取一对数x和y，分别检验x的平方和y的平方之和小于等于1是否成立。设n次试验中等式成立的共有m次，令4mn。但这种方法很难得到的较好的近似值。编写程序：m=0;for i=1:100000x=rand;y=rand;if x2+y2<=1;m=m+1;else end end4*m/100000得3.136000000000000n=100

11、000时，3.1399200000000003. 数值积分方法 半径为1的圆的面积是。以圆心为原点建立直角坐标系，则圆在第一象限的扇形是由与轴，轴所围成的图形，扇形的面积。只要求出扇形的面积，就可得出的值。而扇形面积可近似等于定积分的值。 对于定积分的值，可以看做成曲线与轴，所围的曲边梯形的面积。把分成等分，既得个点，组成个小区间，每一个小区间与轴，所围成的图形是一个小曲边梯形。而梯形的面积计算公式是，对于第个小曲边梯形有上底为，下底为。所有小梯形的高都为。所以第个小曲边梯形的面积为。曲边梯形的总面积即定积分的值就是所有小梯形的面积总和。 为了避免根号，我们也可以利用积分得出的值。我们可以利用

12、对求曲边梯形的面积来得出定积分的值，从而得出的值。设分点x1,x2,xn-1将积分区间0,1分成n等分。所有的曲边梯形的宽度都是h=1/n。记yi=f(xi).则第i个曲边梯形的面积A近似地等于梯形面积，即：A=(y(i-1)+yi)h/2。将所有这些梯形面积加起来就得到：A2/n2(y1+y2+yn-1)+y0+yn编写程序：n=10000;x=linspace(0,1,n+1);y=1./(1+x.2);h=1/n;ans=vpa(4*h*trapz(y),11)得3.1415926519n=100000时，编写程序：n=100000;x=linspace(0,1,n+1);y=1./(1

13、+x.2);h=1/n;ans=vpa(4*h*trapz(y),11)得3.1415926536也可利用积分公式编写程序：n=10000;x=linspace(0,1,n+1);y=sqrt(1-x.2);h=1/n;ans=vpa(4*h*trapz(y),11)得3.1415914776n=100000时，得3.1415926164结果分析对于分析方法来说，一般而言逼近速度太慢，运算庞大，对速度造成了很大影响。相对来说麦琴和泰勒级数法逼近的速度大大增加，得到的近似值精度较高。对于概率方法来说，随机性很大，同一个实验次数，得出的并不相同，有时差别还会很大，所以这种方法很难得到较好的近似值。

14、 对于数值积分法来说，计算量较大，运算慢，不如分析方法中的麦琴和泰勒级数法好，泰勒级数法的精确度较高。作业二：正态分布随机数的产生方法一：平均值为a，标准偏差为b，长度为N的正态分布MATLAB程序如下：m=input('请输入平均值：');n=input('请输入标准差：');t=input('请输入数据长度：');%产生正态分布的随机数for i=1:t a=rand; b=rand; X(i)=sqrt(-2)*log(a)*cos(2*pi*b); Y=X*n+m;enddisp(Y);%求平均值和标准差M=mean(Y);N=std(Y

15、);disp(M);disp(N);%将数据写入文本文件fid=fopen('noise.dat','w');Z=Y;fprintf(fid,'%ft',Z);fclose(fid);%绘图histfit(Y);xlabel('随机数');ylabel('出现的次数');%检验h=lillietest(Y); %若结果h为1，则说明零假设不成立，拒绝零假设；否则， 结果为0，零假设成立，即原分布为正态分布。disp(h); 当N=1000时，Elapsed time is 0.312030 seconds.运行时间

16、在0.3s左右。当N=10000时，Elapsed time is 0.581916 seconds.运行时间在0.6s左右。 当均值a=4，标准偏差b=1，长度N=10000时，得到的图像如下： 实际得到的平均值=4.0069，标准偏差=0.9999，与设定值之间的误差很小，直方图与正态分布曲线基本吻合。需注意这个问题中，我使用的函数为rand随机数发生器，其实matlab中的随机函数并不是真正意义上的随机函数，而是按照一定的递推规则产生的伪随机数。但是在一定的可信度范围内，可以认为是真正的随机数。方法二：利用中心极限定理方法的实现与分析利用中心极限定理来生成随机数的函数(C+语言)编写如下

17、：const int N = 200;double getRand()double s = 0;for (int i = 0; i != N; +i)s += double(rand() % 1000) / 1000;return s; 函数生成的随机数是N个0,1间服从均匀分布的随机数的和。这里N为200。从而理论上产生的随机数应近似服从，其中n为N，即200，为0.5，为1/12。程序生成了200个随机数，并求出样本均值与样本方差，也即与的最大似然估计：/生成随机数并存储double sum, store200, xi, su = 0, sb = 0, ssb = 0;int cnt =

18、0;sum = 0;for (int i = 0; i != 200; +i)xi = getRand();sum += xi;storei = xi;/得到样本均匀与样本方差su = sum / 200;for (int i = 0; i != 200; +i)sb += (storei - su) * (storei - su);sb /= 200;ssb = sqrt(sb); 此次选取，它们将实轴分成11个互不相交的区间，从而将样本值分成11组。 程序统计了每组中的样本数量。为方便计算，程序还计算出了：int segments12, m = 2;double x1 = 90, x10

19、= 108;memset(segments, 0, sizeof(segments);for (int i = 0; i != 200; +i)if (storei <= x1)+segments0;else if (storei > x10)+segments10;else+segmentsint(storei - x1) / m + 1);cout << 'i' << 't' << "ni" << endl;for (int i = 0; i != 11; +i)cout <

20、;< i + 1 << 't' << segmentsi;if (i < 10)cout << 't' << fixed << setprecision(2) << (90 + i * 2 - su) / ssb;cout << endl;程序的最终运行输出如图21所示。图21 最终输出结果对结果的统计如表21所示。由表21中可见，今并令,则由于，故可认为产生的随机数服从正态分布。表21方法三：Box-Muller方法 Box-Muller算法一般是要得到服从正态分布的随机数，基本思想是先得到服从均匀分布的随机数，再将服从均匀分布的随机数转变为服从正态分布。Box-Muller变换通常由2种形式表示，极坐标形式和基本形式。Box给出了基本形式，Muller给出了两个在区间（0，1 上服从均匀分布的样本，再把它们映射为两个标准正态分布的样本。极坐标形式需要两个来自区间 1，1 的不同样本，并将它们映射为两个正态分布的样本，但不使用正弦或余弦函数。 具体地，将X和Y两个标准正态随机变量用极坐标 R和 表示，由独立性，则可写出X，Y的联