統(tǒng)計計算統(tǒng)計模擬π的估計

統(tǒng)計模擬作業(yè)作業(yè)一：n的估計1.分析方法從十七世紀中葉起，人們開始用更先進的分析方法來求n的近似值，其中應用的主要工具是收斂的無窮乘積和無窮級數(shù)，下面列舉一些用此類方法求n近似值的實例。(1)由公式仝二i-仝(1)由公式仝二i-仝(—1)kk=o12k+1推出：=4￡(—1)kk=o12k+1編寫程序:symskx二symsum((T)八k/(2*k+1),k,0,10)y=4*x得出當k=10時，n~3.232315809405593編寫程序symskx二symsum((T)八k/(2*k+1),k,0,20)y=4*x得出當k=20時，n~3.189184782277595依次，加大k的值K=50,n~3.161198612987050K=100,n~3.151493401070990K=200,n~3.146567747182986Kn103.232315809405593503.1611986129870501003.1514934010709902003.146567747182986沃里斯(Wallis)方法在積分學中我們經(jīng)常會遇到如下的沃利斯(Wallis)公式。沃利斯(Wallis)公式揭示了n與整數(shù)之間的一種很不尋常的關系。但在實際學習中很少注意到沃利斯(Wallis)公式,更不會關注它的應用。實際上,沃利斯(Wallis)公式有許多作用,經(jīng)常有以下幾方面的應用。1應用于極限計算中。由于沃利斯(Wallis)公式與極限有關，所以有些極限的計算可以通過沃利斯(Wallis)公式很容易計算出來。2.應用于積分計算中。對于一些用積分法不易求得原函數(shù)的積分，而用沃利斯(Wallis)公式卻很容易解決問題。編寫程序：formatlongx=1;fork=1:10x=x*(2*k/(2*k-1)*2*k/(2*k+1));endy=2*x得k=10時，n~3.067703806643498增加k的值K=20,n~3.103516961539230K=50,n~3.126078900215409K=100,n~3.133787490628159K=10000,n~3.141514118681864K=1000000,n~3.141591868191880kn103.067703806643498203.103516961539230503.1260789002154091003.133787490628159100003.14151411868186410000003.141591868191880（3）蒙特卡羅法取一正方形A,以A的一個頂點為圓心，A的邊長為半徑畫圓，取四分之一圓（正方形內(nèi)的四分之一圓）為扇形B。已知A的面積，只要求出B

的面積與A的面積之比的面積與A的面積之比k=SBA，就能得出S，B再由B的面積為圓面積的四分之一，利用公式S=,R2即可求出“的值。因此，我們的目的就是要找出k圓的值?？梢园袮和B看成是由無限多個點組成，而B內(nèi)的所有點都在A內(nèi)。隨機產(chǎn)生n個點，若落在B內(nèi)的有m個點（假定A的邊長為1,以扇形圓心為坐標系原點。則只要使隨機產(chǎn)生橫縱坐標x、y滿足x2+y2<1的點，就是落在b內(nèi)的點），則可近似得出k的值，即k二m，由此就可以求出“的值。n編寫程序：i=1;m=0;n=1000;fori=1:na=rand(1,2);ifa(1廠2+a(2廠2<=1m=m+1;endendp=vpa(4*m/n,30)程序運行結(jié)果：p=3.14000000000000000000000000000泰勒級數(shù)法函數(shù)arctanx的泰勒級數(shù)展開式為x3 x5 x2k—1arctanx=x— + +(―1)k—1 +…3 5 2k—1將x=1代入上式有兀 1 1 1 1 / 1、 1=arctanl=1——+——???+(―1)n-1—4 35 2n—1利用這個式子就可以求出兀的值了。編寫程序：i=1;n=1000;s=0;fori=1:ns二s+(-1廠(i-1)/(2*i-1);endp=vpa(4*s,30)程序運行結(jié)果P=3.14059265383979413499559996126當取n的值為10000時，就會花費很長時間，而且精度也不是很高。原因是x=1時，arctan1的展開式收斂太慢。因此就需要找出一個x使得arctanx收斂更快。若取x=丄，則我們只有找出?與上的關系，才能求出兀的值。2 4令 arctan—，P=^—a，根據(jù)公式tan(a+P)=tan^+tan弓有tanP=1，則有2 4 1—tanatanP 3parctan2+arctan|。所以可以用專二耐專+arctan|來計算兀的值。

利用公式L=arctan1=arctan丄+arctan—4 2 3推出n=4(arctan2 3編寫程序：i=1;n=1000;s=0;s1=0;s2=0;fori=1:nsi二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*s,30)程序運行結(jié)果p=3.14159265358979323846264338328TOC\o"1-5"\h\z顯然，級數(shù)收斂越快，取同樣的n值可以得到更高的精度。以同樣的方法，能得出L=4arctan-+arctan丄，程序和上面的一樣。這樣"的近似值可4 5 239以精確到幾百位。6）麥琴公式由上=4arctan-arctan—這個公式由英國天文學教授約翰?這個公式由英國天文學教授約翰?麥琴于推出n=4(4arctan—-arctan」一),5 2391706年發(fā)現(xiàn)。他利用這個公式計算到了100位的圓周率。麥琴公式每計算一項可以得到1.4位的十進制精度。因為它的計算過程中被乘數(shù)和被除數(shù)都不大于長整數(shù)，所以可以很容易地在計算機上編程實現(xiàn)。編寫程序：symsn;= n-1)*(1/5廠(2*n-1)/(2*n-1);f2=(-1)^(n-1)*(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)用下列公式求n值，并了解n的其他公式。編寫程序：digits(10000);N=100000000;p=1;form=1:N;m1=m*m;m2=2*m;n=(4*m1)/((m2-1)*(m2+1));p=p*n;endp=vpa(p*2);改變N的值，分析得到的不同n值的大小與N的關系。PI=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825當N=10*10八4時，n=3.14151411868195662435709891724400222301483154296875當N=10*10八8時，n=3.141592643066262180440162410377524793148040771484375在使用MATLAB的過程中發(fā)現(xiàn)，當N值越大的時候，PI的值精確度越高，然而N大到一定范圍后得出結(jié)果的過程也就變得非常慢，并且在10*10八8后每增加10倍后PI沒有很大的變化，并且運算速度也變得非常慢了。2.概率方法取一個二維數(shù)組（x,y）,取一個充分大的正整數(shù)n,重復n次，每次獨立的從（0,1）中隨機地取一對數(shù)x和y,分別檢驗x的平方和y的平方之和小于等于1是否成立。設n次試驗中等式成立的共有m次,令n~4m/n。但這種方法很難得到n的較好的近似值。編寫程序：m=0;fori=1:100000x=rand;y=rand;ifx八2+y八2<=1;m=m+1;elseendend4*m/100000得n~3.136000000000000n=100000時，n~3.139920000000000數(shù)值半徑為1的圓的面積是兀。以圓心為原點建立直角坐標系，則圓在第—象限的扇形是由yk'G與x軸，y軸所圍成的圖形，扇形的面積s冷。只要求出扇形的面積，就可得出兀的值。而扇形面積可近似等于定積分一X2dx的值。0對于定積分Jbf(x)dx的值，可以看做成曲線與X軸，X二a,x=b所圍的a曲邊梯形的面積S。把［a,b］分成n寺分，既得n+1個點x=a,x,-x,x=b1 n一1組成n個小區(qū)間，每一個小區(qū)間與x軸，x=a,x=b所圍成的圖形是一個小曲邊梯形。而梯形的面積計算公式是(上底+下底)x高一2，對于第i個小曲邊梯形有上底為f(x),下底為/(x)。所有小梯形的高都為h=(b-a)/n。所以i一1 i第i個小曲邊梯形的面積為(f(x)+f(x))xh十2。曲邊梯形的總面積即定積分i一1 iJbf(x)dx的值就是所有小梯形的面積總和。a為了避免根號，我們也可以利用積分j1 -dx=o1+x2 4得出兀的值。我們可以利用對求曲邊梯形的面積來得出定積分j1丄dx的值，從而得出兀o1+x2的值。設分點x,x,…X將積分區(qū)間[0,1]分成n等分。1 2 n-1所有的曲邊梯形的寬度都是h=1/n。記yi二f(xi).則第i個曲邊梯形的面積A近似地等于梯形面積，即:A=(y(i-1)+yi)h/2。將所有這些梯形面積加起來就得到：A~2/n[2(y+y+…y)+y+y]1 2 n-1 0n編寫程序：n=10000;x=linspace(0,1,n+1);y=1./(1+x.八2);h=1/n;ans=vpa(4*h*trapz(y),11)得n~3.1415926519n=100000時，編寫程序：n=100000;x=linspace(0,1,n+1);y=1./(1+x.八2);h=1/n;ans=vpa(4*h*trapz(y),11)3.1415926536也可利用積分公式“二4卩訂二2dx0編寫程序：n=10000;x=linspace(0,1,n+1);y二sqrt(1-x.八2);h=1/n;ans=vpa(4*h*trapz(y),11)得n~3.1415914776n=100000時，得n~3.1415926164結(jié)果分析對于分析方法來說，一般而言逼近速度太慢，運算龐大，對速度造成了很大影響。相對來說麥琴和泰勒級數(shù)法逼近的速度大大增加，得到的近似值精度較高。對于概率方法來說，隨機性很大，同一個實驗次數(shù)，得出的n并不相同，有時差別還會很大，所以這種方法很難得到n較好的近似值。對于數(shù)值積分法來說，計算量較大，運算慢，不如分析方法中的麥琴和泰勒級數(shù)法好，泰勒級數(shù)法的精確度較高。作業(yè)二：正態(tài)分布隨機數(shù)的產(chǎn)生方法一：平均值為a,標準偏差為b,長度為N的正態(tài)分布MATLAB程序如下：m=input('請輸入平均值：’)；n=input('請輸入標準差:');t=input('請輸入數(shù)據(jù)長度：’)；％產(chǎn)生正態(tài)分布的隨機數(shù)fori=1:ta=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);disp(M);disp(N);%將數(shù)據(jù)寫入文本文件fid=fopen('noise.dat','w');Z=Y;fprintf(fid,'%f\t',Z);fclose(fid);%繪圖histfit(Y);xlabel('隨機數(shù)');ylabel('出現(xiàn)的次數(shù)');%檢驗h=lillietest(Y);%若結(jié)果h為1，則說明零假設不成立，拒絕零假設；否則，結(jié)果為0,零假設成立，即原分布為正態(tài)分布。disp(h);當N=1000時，Elapsedtimeis0.312030seconds■運行時間在0.3s左右。當N=10000時，Elapsedtimeis0.581916seconds.運行時間在0.6s左右。當均值a=4,標準偏差b=1，長度N=10000時，得到的圖像如下：蘇長蘇長S皋刃實際得到的平均值=4.0069，標準偏差=0.9999,與設定值之間的誤差很小，直方圖與正態(tài)分布曲線基本吻合。需注意這個問題中，我使用的函數(shù)為rand隨機數(shù)發(fā)生器，其實matlab中的隨機函數(shù)并不是真正意義上的隨機函數(shù)，而是按照一定的遞推規(guī)則產(chǎn)生的偽隨機數(shù)。但是在一定的可信度范圍內(nèi)，可以認為是真正的隨機數(shù)。方法二：利用中心極限定理方法的實現(xiàn)與分析利用中心極限定理來生成隨機數(shù)的函數(shù)(C++語言)編寫如下:constintN=200;doublegetRand(){doubles=0;for(inti=0;i!=N;++i)s+=double(rand()%1000)/1000;returns;}函數(shù)生成的隨機數(shù)是N個[0,1]間服從均勻分布的隨機數(shù)的和。這里N為200。從而理論上產(chǎn)生的隨機數(shù)應近似服從n(ny,nc2)，其中n為N，即200，卩為0.5,°2為1/12。程序生成了200個隨機數(shù)，并求出樣本均值與樣本方差，也即卩與°2的最大似然估計：//生成隨機數(shù)并存儲doublesum,store[200],xi,su=0,sb=0,ssb=0;intcnt=0;sum=0;for(inti=0;i!=200;++i){xi=getRand();sum+=xi;store[i]=xi;}//得到樣本均勻與樣本方差su=sum/200;for(inti=0;i!=200;++i)sb+=(store[i]-su)*(store[i]-su);sb/=200;ssb=sqrt(sb);此次選取y二90,y二92,y二94,y二96,y二9&y二100,y二102,y二104,1 2 3 4 5 6 7 810y二106,y二108，它們將實軸分成11個互不相交的區(qū)間，從而將樣本值分成91011組。程序統(tǒng)計了每組中的樣本數(shù)量。為方便計算，程序還計算出了(yi?7):6intsegments[12],m=2;doublex1=90,x10=108;memset(segments,0,sizeof(segments));for(inti=0;i!=200;++i){if(store[i]<=x1)++segments[0];elseif(store[i]>x10)++segments[10];else++segments[int((store[i]-x1)/m+1)];}cout<<'i'<<'\t'<<"ni"<<endl;for(inti=0;i!=11;++i){cout<<i+1<<'\t'<<segments[i];if(i<10)cout<<'\t'<<fixed<<setprecision(2)<<(90+i*2-su)/ssb;cout<<endl;}程序的最終運行輸出如圖2－1所示。

彳羊本均值：丫9?2773片羊本方差：16?8202k：lln-2ini.11-2-2625-1-77315-1-29424-0-80529-0-3163&0.1S73?0.6682?1-159171-641032-13114圖2—1最終輸出結(jié)果對結(jié)果的統(tǒng)計如表2—1所示。由表2—1中可見咒2二13.380，今k二11,m二2,并令a=0.05,則%2(k-m-a1)=%2(8)=15.507.由于13.380<15.507，故可認為產(chǎn)生0.05的隨機數(shù)服從正態(tài)分布。i(y,y]i-1 inipi(n-200p)2?200p?i1(-8,90]20.01190.0612(90,92]20.02652.0553(92,94]100.06010.0344(94,96]120.11345.0295(96,98]370.16640.4166(9&100]460.20690.5167(100,102]420.17401.4898(102,104]230.12950.3259(104,106]150.07460.00010(106,108]100.03391.52911(108,+刈10.01661.62113.380表21方法三：Box-Muller方法Box-Muller算法一般是要得到服從正態(tài)分布的隨機數(shù)，基本思想是先得到服從均勻分布的隨機數(shù)，再將服從均勻分布的隨機數(shù)轉(zhuǎn)變?yōu)榉恼龖B(tài)分布。Box-Muller變換通常由2種形式表示，極坐標形式和基本形式。Box給出了基本形式，Muller給出了兩個在區(qū)間(0,1］上服從均勻分布的樣本，再把它們映射為兩個標準正態(tài)分布的樣本。極坐標形式需要兩個來自區(qū)間［1,1］的不同樣本，并將它們映射為兩個正態(tài)分布的樣本，但不使用正弦或余弦函數(shù)。具體地，將x和丫兩個標準正態(tài)隨機變量用極坐標r和e表示，由獨立性，則可寫出X,Y的聯(lián)合密度，從而得到和0的聯(lián)合密度。觀察此聯(lián)合密度的形式，得TT和㈢對應的特殊分布。下面則是生成標準正態(tài)隨機變量的Box-Muller算法：步驟一：生成隨機數(shù)U1U2.步驟二：R2=—21nU(人」是均值為2的指數(shù)隨機變量)?二2兀U2(。是(0,12:-)上均勻隨機變量)步驟三：令X=Rcos0=—21nUcos(2兀U) Y=Rsin?=「一21nUsin(2兀U)1212步驟三中兩式的變換稱為Box-Muller變換使用BoxMuller方法得到隨機數(shù)的函數(shù)如下：doublegetRand(){doubleu1=double(rand()%10000)/10000,u2=double(rand()%10000)/10000,r;r=20+5*sq