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

统计模拟作业班级学号姓名
作业一：π的估计
1. 分析方法
从十七世纪中叶起，人们开始用更先进的分析方法来求∏的近似值，其中应用的主要工具是收敛的无穷乘积和无穷级数，下面列举一些用此类方法求∏近似值的实例。

(1)由公式121
)1(5131140
+-=-+-=∑∞
=k k k π
推出=4121
)1(0+-∑∞
=k k k
编写程序： syms k
x=symsum((-1)^k/(2*k+1),k,0,10) y=4*x
得出当k=10时，π≈3.232315809405593 编写程序 syms k
x=symsum((-1)^k/(2*k+1),k,0,20) y=4*x
得出当k=20时，π≈3.189184782277595 依次，加大k 的值
K=50，π≈3.161198612987050 K=100，π≈3.151493401070990
K=200，π≈3.146567747182986
…
K π
10 3.232315809405593
50 3.161198612987050
100 3.151493401070990
200 3.146567747182986
（2）沃里斯(Wallis)方法
在积分学中我们经常会遇到如下的沃利斯(Wallis)公式。

沃利斯(Wallis)公式揭示了π与整数之间的一种很不寻常的关系。

但在实际学习中很少注意到沃利斯(Wallis)公式,更不会关注它的应用。

实际上,沃利斯(Wallis)公式有许多作用,经常有以下几方面的应用。

1应用于极限计算中。

由于沃利斯(Wallis)公式与极限有关,所以有些极限的计算可以通过沃利斯(Wallis)公式很容易计算出来。

2.应用于积分计算中。

对于一些用积分法不易求得原函数的积分，而用沃利斯(Wallis)公式却很容易解决问题。

编写程序：
format long
x=1;
for k=1:10
x=x*(2*k/(2*k-1)*2*k/(2*k+1));
end
y=2*x
得k=10时，π≈3.067703806643498
增加k的值
K=20，π≈3.103516961539230
K=50，π≈3.126078*********
K=100，π≈3.133787*********
K=10000，π≈3.141514118681864
K=1000000，π≈3.141591868191880
…
k π
10 3.067703806643498
20 3.103516961539230
50 3.126078*********
100 3.133787*********
10000 3.141514118681864
1000000 3.141591868191880
（3）蒙特卡罗法
取一正方形A，以A的一个顶点为圆心，A的边长为半径画圆，取四分之一圆（正方形内的四分之一圆）为扇形B。

已知A的面积，只要求出B
的面积与A 的面积之比B
A
S k S =
，就能得出B S ，再由B 的面积为圆面积的四分之一，利用公式2=S R π圆即可求出π的值。

因此，我们的目的就是要找出k 的值。

可以把A 和B 看成是由无限多个点组成，而B 内的所有点都在A 内。

随机产生n 个点，若落在B 内的有m 个点（假定A 的边长为1，以扇形圆心为坐标系原点。

则只要使随机产生横纵坐标x 、y 满足221x y +≤的点，就是落在B 内的点），则可近似得出k 的值，即m
k n
=，由此就可以求出π的值。

编写程序： 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; end end
p=vpa(4*m/n,30) 程序运行结果： p =
3.14000000000000000000000000000
（4）泰勒级数法
函数arctan x 的泰勒级数展开式为：
3521
1arctan (1)3521
k k x x x x x k --=-+-⋅⋅⋅+-+⋅⋅⋅- 将1x =代入上式有
1111
arctan11(1)43521
n n π
-==-+-⋅⋅⋅+--. 利用这个式子就可以求出π的值了。

编写程序： i=1;n=1000;s=0; for i=1:n
s=s+(-1)^(i-1)/(2*i-1); end
p=vpa(4*s,30) 程序运行结果： p =
3.14059265383979413499559996126
当取n 的值为10000时，就会花费很长时间，而且精度也不是很高。

原因是1x =时，arctan1的展开式收敛太慢。

因此就需要找出一个x 使得arctan x 收敛更快。

若取1
2x =，则我们只有找出α与4
π的关系，才能求出π的值。

令1arctan 2
α=，4
π
βα=
-，根据公式tan tan tan()1tan tan αβ
αβαβ
++=
-有1tan 3β=，则有
11arctan
arctan 4
23π
=+。

所以可以用11
arctan arctan 423
π=+来计算π的值。

（5）利用公式
3
1arctan 21arctan
1arctan 4
+==π
推出π=4(3
1
arctan 21arctan +) 编写程序：
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); end s=s1+s2; p=vpa(4*s,30) 程序运行结果： p =
3.14159265358979323846264338328
显然，级数收敛越快，取同样的n 值可以得到更高的精度。

以同样的方法，能得出11
4arctan arctan
4
5239
π
=+，程序和上面的一样。

这样π的近似值可以精确到几百位。

（6）麦琴公式
由
239
1arctan 51arctan
44-=π
推出π=4(239
1
arctan 51arctan 4-)，这个公式由英国天文学教授约翰·麦琴于
1706年发现。

他利用这个公式计算到了100位的圆周率。

麦琴公式每计算一项可以得到1.4位的十进制精度。

因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。

编写程序：
syms n;
f1=(-1)^(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.14159265358979323846264338327950288419713045104626857897220325 5663716036677133432949011735665451127
（7）用下列公式求π值，并了解π的其他公式。

编写程序：
digits(10000);
N=100000000;p=1;
for m=1:N;m1=m*m;m2=2*m;
n=(4*m1)/((m2-1)*(m2+1));
p=p*n;
end
p=vpa(p*2);
改变N的值，分析得到的不同π值的大小与N的关系。

PI=3. 1415926535897932384626433832795028841971693993751058209749 44592307816406286208998628034825
当N=10*10^4时，
π=3.14151411868195662435709891724400222301483154296875
当N=10*10^8时，
π=3.141592643066262180440162410377524793148040771484375 在使用MATLAB的过程中发现，当N值越大的时候，PI的值精确度越高，然而N大到一定范围后得出结果的过程也就变得非常慢，并且在10*10^8后每增加10倍后PI没有很大的变化，并且运算速度也变得非常慢了。

2.概率方法
取一个二维数组 (x,y) ,取一个充分大的正整数n,重复n次，每次独立的从（0,1）中随机地取一对数x和y，分别检验x的平方和y的平方之和小于等于1是否成立。

设n次试验中等式成立的共有m次，令∏≈4m／n。

但这种方法很难得到∏的较好的近似值。

编写程序：
m=0;
for i=1:100000
x=rand;
y=rand;
if x^2+y^2<=1;
m=m+1;
else end
end 4*m/100000
得π≈3.136000000000000 n=100000时，
π≈3.139920000000000
3.数值积分方法
半径为1的圆的面积是π。

以圆心为原点建立直角坐标系，则圆在第
一象限的扇形是由2
1y x =-与x 轴，y 轴所围成的图形，扇形的面积4s π
=。

只要求出扇形的面积，就可得出π的值。

而扇形面积可近似等于定积分
1
20
1x dx -⎰
的值。

对于定积分()b
a f x dx ⎰的值，可以看做成曲线与x 轴，x a =，x
b =所围的曲边梯形的面积S 。

把[],a b 分成n 等分，既得1n +个点x a =，1x ，⋅⋅⋅1n x -，x b =组成n 个小区间，每一个小区间与x 轴，x a =，x b =所围成的图形是一个小曲边梯形。

而梯形的面积计算公式是(+)2⨯÷上底下底高，对于第i 个小曲边梯形有上底为1()i f x -，下底为()i f x 。

所有小梯形的高都为()/h b a n =-。

所以第i 个小曲边梯形的面积为1(()+())2i i f x f x h -⨯÷。

曲边梯形的总面积即定积分
()b
a
f x dx ⎰
的值就是所有小梯形的面积总和。

为了避免根号，我们也可以利用积分
1
2
0114
dx x π
=+⎰ 得出π的值。

我们可以利用对求曲边梯形的面积来得出定积分1
2
1
1dx x
+⎰的值，从而得出π
的值。

设分点x1,x2,…x n-1将积分区间[0,1]分成n等分。

所有的曲边梯形的宽度都是h=1/n。

记yi=
f(xi).则第i个曲边梯形的面积A近似地等于梯形面积，即：A=(y(i-1)+yi)h/2。

将所有这些梯形面积加起来就得到：
A≈2/n[2(y1+y2+…y n-1)+y0+y n]
编写程序：
n=10000;
x=linspace(0,1,n+1);
y=1./(1+x.^2);
h=1/n;
ans=vpa(4*h*trapz(y),11)
得π≈3.1415926519
n=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
也可利用积分公式⎰-
=1
02
1
4dx
x
π
编写程序：
n=10000;
x=linspace(0,1,n+1);
y=sqrt(1-x.^2);
h=1/n;
ans=vpa(4*h*trapz(y),11)
得π≈3.1415914776
n=100000时，
得π≈3.1415926164
结果分析
对于分析方法来说，一般而言逼近速度太慢，运算庞大，对速度造成了很大影响。

相对来说麦琴和泰勒级数法逼近的速度大大增加，得到的近似值精度较高。

对于概率方法来说，随机性很大，同一个实验次数，得出的π并不相同，有时差别还会很大，所以这种方法很难得到π较好的近似值。

对于数值积分法来说，计算量较大，运算慢，不如分析方法中的麦琴和泰勒级数法好，泰勒级数法的精确度较高。

作业二：正态分布随机数的产生
方法一：
平均值为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;
end
disp(Y);%求平均值和标准差
M=mean(Y);
N=std(Y);
disp(M);
disp(N);%将数据写入文本文件
fid=fopen('noise.dat','w');
Z=Y;
fprintf(fid,'%f\t',Z);
fclose(fid);%绘图
histfit(Y);
xlabel('随机数');
ylabel('出现的次数');%检验
h=lillietest(Y); %若结果h为1，则说明零假设不成立，拒绝零假设；否则，结果为0，零假设成立，即原分布为正态分布。

disp(h);
当N=1000时，Elapsed time is 0.312030 seconds.运行时间在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++语言)编写如下：
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。

从而理论上产生的随机数应近似服从2
μσ，其中n为N，即
N n n
(,)
200，μ为0.5，2σ为1/12。

程序生成了200个随机数，并求出样本均值与样本方差，也即μ与2σ的最大似然估计：
//生成随机数并存储
double sum, store[200], xi, su = 0, sb = 0, ssb = 0;
int cnt = 0; sum = 0;
for (int i = 0; i != 200; ++i) {
xi = getRand(); sum += xi; store[i] = xi; }
//得到样本均匀与样本方差 su = sum / 200;
for (int i = 0; i != 200; ++i)
sb += (store[i] - su) * (store[i] - su); sb /= 200; ssb = sqrt(sb);
此次选取1234567890,92,94,96,98,100,102,104,y y y y y y y y ======== 910106,108y y ==，它们将实轴分成11个互不相交的区间，从而将样本值分成11组。

程序统计了每组中的样本数量。

为方便计算，程序还计算出了ˆ()ˆi y μ
σ
-： int segments[12], m = 2; double x1 = 90, x10 = 108;
memset(segments, 0, sizeof (segments));
for (int i = 0; i != 200; ++i)
{
if (store[i] <= x1)
++segments[0];
else if (store[i] > x10)
++segments[10];
else
++segments[int((store[i] - x1) / m + 1)];
}
cout << 'i' << '\t' << "ni" << endl;
for (int i = 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所示。

图2－1 最终输出结果
对结果的统计如表2－1所示。

由表2－1中可见213.380χ=，今11,2,k m ==并
令0.05α=,则()()22
0.051815.507.k m αχχ--==由于13.38015.507<，故可认为产生
的随机数服从正态分布。

表2－1
方法三：Box-Muller 方法
Box-Muller 算法一般是要得到服从正态分布的随机数，基本思想是先
i
1(,]i i y y -
i n
ˆi p
()2
ˆ200ˆ200i i i n p p
-
123
4567891011
(,90](90,92](92,94](94,96](96,98](98,100](100,102](102,104](104,106](106,108](108,]
-∞+∞ 2210123746422315101
0.01190.02650.06010.11340.16640.20690.17400.12950.07460.03390.0166
0.0612.0550.0345.0290.416
0.5161.4890.3250.0001.5291.621 ∑
13.380
得到服从均匀分布的随机数，再将服从均匀分布的随机数转变为服从正态分布。

Box-Muller 变换通常由2种形式表示，极坐标形式和基本形式。

Box 给出了基本形式，Muller 给出了两个在区间（0，1 ]上服从均匀分布的样本，再把它们映射为两个标准正态分布的样本。

极坐标形式需要两个来自区间 [− 1，1 ]的不同样本，并将它们映射为两个正态分布的样本，但不使用正弦或余弦函数。

具体地，将X 和Y 两个标准正态随机变量用极坐标 R 和 θ表示，由独立性，则可写出X ，Y 的联合密度，从而得到和的联合密度。

观
察此联合密度的形式，得
和对应的特殊分布。

下面则是生成标准正态随机变量的Box-Muller 算法： 步骤一：生成随机数U1 U2. 步骤二：12
l 2nU R -=（是均值为2的指数随机变量），U 22π=Θ（是（0,
）上均匀随机变量）
步骤三：令
)2cos(ln 2cos R X 21U U π-=Θ= )U 2sin(lnU 2sin R Y 21π-=Θ=
步骤三中两式的变换称为Box-Muller 变换 使用Box Muller 方法得到随机数的函数如下： double getRand() {
double u1 = double (rand() % 10000) / 10000, u2 = double (rand() % 10000) / 10000, r;
r = 20 + 5 * sqrt(-2.0 * (log(u1) / log(e))) * cos(2 * pi * u2); return r; }
用此函数得到的随机数样本理论上服从(20,25)N 。

所实现的程序产生了500个随机变量的样本其他与利用中心极限定理的实现基本相同。

程序得到的最终结果如图2－2所示。

图2-2
对结果的统计如表2-2 所示。

表2－2
i
1(,]i i y y -
i n ˆi p
()2
ˆ500ˆ500i i i n p p
-
２０
1234567891011
(,10](10,12](12,14](14,16](16,18](18,20](20,22](22,24](24,26](26,28](28,]
-∞+∞ 1711354758649261553030
0.02500.03320.06080.09580.12980.15140.15200.13140.09860.07240.0606
1.621.890.690.010.731.813.360.330.651.060.00
∑
12.19
由表
2－2
可见212.19χ=，今11,2,k m ==并令0.0α=,则
()()22
0.051815.507.k m αχχ--==由于12.1915.507<，
故可认为产生的随机数服从正态分布。

结果分析
利用中心极限定理的方法虽然可以得到服从正态分布的随机数样本，其思想也较为简单，容易想到。

但是这种方法每次都要先产生若干个服从均匀分布的随机数样本并求它们的和，因而算法的时间复杂度高。

Box Muller 方法的推导过程较为复杂。

但得到的结果却是很令人满意的。

只要用两个相互独立的均匀分布就能得到正态分布。

而且产生随机数的时间复杂度比利用中心极限定理的方法要低很多。

因而若要产生服从正态分布的随机数样例，则Box Muller 方法是一个很不错的选择。