哈工大数学实验大作业 Galton钉板问题

1、Galton钉板问题一、实验目的与要求1.复习概率论中随机变量、概率分布、二项分布、均值和分布函数等概念。2.理解Galton钉板实验中小球落入格子所服从的规律。3.了解Matlab软件中进行动画演示的命令。4.掌握Matlab软件中进行随机模拟的方法。二、问题描述所有现象的“因”和“果”，即“条件”和“结果”之间在客观上都存在着一定的规律，这种规律通常可以分成两类：一是确定性的规律，另一类是非确定性规律。对于确定性的系统，当已知条件是充分时，那么实验的结果也是确定的，即在每一次试验进行以前，可以预见试验产生的结果。但若条件不充分时，就无法预测试验的结果，这就产生了“因果律的缺失”的随机现象。

2、随机现象在实践中是大量遇到的，如掷骰子。虽然无法由“因”预测“果”，但是当进行大量重复试验时，因果之间仍会呈现一种统计规律。概率方法建立在“重复试验”的基础上，统计规律只有在大量重复后才会呈现出来，诸如随机变量。分布、均值、方差等概念无一不体现了重复的思想。以下围绕着Galton钉板模型来讨论。Galton钉板试验是由英国生物统计学家Galton设计的。在一板上有n排钉子，图0所示的是n4的情况。图中10个圆点表示10颗钉子，在钉子的下方有5个格子，分别编号为0，1，2，3，4。自Galton钉板的上方扔进一个小球任其自由下落，在下落的过程中当小球碰到钉子时，从左边落下与从右边落下的机会相等、

3、碰到下一排钉子时又是如此。最后落入底板中的某一个格子。图0 Galton钉板模型（n4）问题：1、一个一个从顶部放入k个小球，底槽中各格的理论频数应为多少？2、（1）编写一演示程序（用MATLAB语言），其界面与完成的基本功能如书中P121.图8.9所示； （2）利用自由度为4的2分布临界值表及统计量（其中vi为04号格子中第i个实验频数）。检验假设H0：所实验的分布服从B(4,p)3、用类似方法模拟Poisson分布与几何分布。三、问题分析问题1、小球自上方落下，经过n个钉子。每经过一个钉子时只有两种可能结果：向右或向左。这是一个具有两个结果（成功和失败）的随机试验E，将向右视为成功，成功的

4、概率为p，向左为失败，失败的概率为q=1-p。小球碰到一个钉子下落一格，相当于进行了一次试验E。小球自顶端落下，碰到n个钉子，最终落在某个格子的过程，恰好相当于将试验E重复了n次，因此一次投球过程就是一个n重贝努利试验（将仅有两个相互排斥结果的试验E独立重复n次，构成了n重贝努利试验En）。n重贝努利试验的成功次数X正好是小球向右移动的次数，它是一个随机变量。根据概率论的结果有XB(n,p)。对于一个随机变量，我们首先要弄清楚它的取值范围，X的取值范围为0,1,2,n，这是什么意思呢？在Galton钉板模型中X0表示小球向右移动的次数，也就是小球一直向左移动，所以它恰好要落在编号为0的格子里；

5、同理X1表示小球恰好要落在编号为1的格子里，依次类推，这就是说，X是小球最终落进的格子编号数，当然它也对应为小球向右移动的次数。二项随机变量的分布列为： 则在Galton钉板实验中，p0.5，底槽中各格的频数应为k*pi。问题2、演示程序已经完成，具体程序会在下文“实验过程”中体现，以下是演示结果如图1：图1 Galton钉板问题演示结果可以从界面看到，我们的演示程序能方便地对k和p的值进行修改再演示，同时还能统计落入每个格子的小球数。在演示界面下方，有一列是用来检验2上侧分位数是否满足等式P22(n)=我们在程序中设置了=(0.005,0.01,0.025,0.05,0.10,0.25,0.

6、75,0.90,0.95,0.975,0.99)等12个值进行检验，发现只有前6个数可以通过检验。同时，在演示界面右下方有一个“HELP”按键，点击它会打开一个word文档，是本程序的使用说明。如图2图2 HELP文档问题3、Poisson分布与几何分布演示结果如图3和图4：图3 Galton模拟Poisson分布演示结果图4 Galton模拟几何分布演示结果鉴于模拟的难度，我们在演示程序中设定了参数不可更改，尽管可能导致结果有偏差，但我们经过多次演示，发现偏差极小。四、背景知识介绍1、随机变量随机变量是随机试验结果的函数，其特点是在试验前，并不能预知这个函数将取何值，这要凭机会，就是“随机”

7、的意思。一旦试验后，取值就确定了。例如，我在3月31日买了一张奖券，到6月30日开奖。当我买这张奖券时，有人可以对我说：“你中奖的金额是个随机的变量，其值在6月30日抽奖试验做过之后才能确定。”明白了这一点就不难举出许多随机变量的例子。例如，某出租车公司的电话订车中心，一天内接到的订车电话的次数；某射手对一活动靶进行射击，到击中目标为止，所进行的射击次数；从一批灯泡中，任取一只，测定这只灯泡的寿命，等等，这些都是随机变量。2、n重Bernoulli试验当依照一定的质量标准，从大批产品中抽出一件进行产品质量合格性检查时，得出的结果可以是二者之一：“这是件不合格产品”或“这是件合格产品”。如果将这

8、样的抽检一件产品看作是进行一次试验，则试验的结果可以是发生A（这是件不合格产品）或Error! Reference source not found.Error! Reference source not found.（A不发生，即产品质量合格）。称这种只有两个可能结果A（称“成功”）或（称“失败”）的试验为Bernoulli试验。有很多试验，其可能的结果不止两个，但由于人们常只对试验是否发生某一种特定结果感兴趣，因而可将之归结为Bernoulli试验、例如，明天的天气可以有多种情况，但若只关心明天是否下雨，则观察明天的天气（作为一次试验），其结果只有两个：“下雨”或“不下雨”，因而可被看作是

9、个Bernoulli试验。实际上常要考察独立重复进行Bernoulli试验的序列，并将这一独立重复试验序列作为单独一个复合试验来对待。这里，所谓独立重复进行Bernoulli试验的意思是，这个序列中的每一试验的结果都只能发生A或，且发生A的概率一样，是某个值p=P(A)。当然，发生的概率也就是一样地是q=P()，并且每一试验发生的结果不会影响其他试验出现的结果。独立重复试验序列最重要的特性是序列由独立重复进行n次Bernoulli试验组成，简称为n重Bernoulli试验。在很多问题中可以用上n重Bernoulli试验模型。例如，若学校的电话总机设有99个分机，已知每号分机平均每小时有3分钟要

10、使用外线，在考虑该总机应设置多少条外线合适的问题时，可归结为n重Bernoulli试验的问题。在任一时刻考察一部分机是否占用外线时，其可能结果只有两个：“占用”（发生A）、“不占用”（发生）。而且据已知数据有p=P(A)=3/60=0.05，所以这是一个p=0.05的Bernoulli试验。由于各分机是否在占用外线可合理地认为是相互独立的，因而这个问题可看成涉及了一个p=0.05的99重Bernoulli试验。再比如，已知某疾病的发病率为0.001，当卫生部门要对一个拥有5000名员工的单位估计此种疾病的发病情况时，需用p=0.001的n重Bernoulli试验模型，这里n5000。3、二项分

11、布在“成功”概率是p，即p=P(A)的n重Bernoulli试验中，事件A出现的次数是二项分布随机变量，其可能的取值为：0，1，n有分布律这个值也被记作b(k;n,p)，服从参数为n，p的二项分布，也记作B(n,p)。4、离散型随机变量的数学期望设随机变量具有概率分布列则当时，称为随机变量的数学期望或均值，记作数学期望表征的是随机变量取值的“平均值”。五、实验过程1、Matlab命令简介命令功能rand(m,n)产生mn个(0,1)区间中的随机数，并将这些随机数存于一个mn矩阵中。每次调用rand(m,n)的结果都会不同。rand(seed,s)如果想保持结果一致，可与rand(seed,s)

12、配合使用，s为一正整数。Moviein(n)创建动画矩阵，制作动画矩阵数据Getframe拷贝动画矩阵Movie(mat,m)播放动画矩阵m次例如：输入命令：rand(seed,1), u=rand(1,6)得到结果：0.5129 0.4605 0.3504 0.0950 0.4337 0.7092而且再次运行该指令时结果保持不变，除非重新设置种子seed的值。例如：输入命令：rand(seed,2),u=rand(1,6) 得到结果：u=0.0258 0.9210 0.7008 0.1901 0.8673 0.4185这样结果才会产生变化。2、动画模拟Galton钉板试验运行观察程序main

13、.fig，屏幕将出现图形窗口，动画模拟扔球过程。如图5是模拟向一个4层Galton钉板扔500次小球的过程的最后的结果。在模拟过程中我们看到，每一个小球落在哪一个格子是无法预测的，但小球逐渐堆积成一种单峰的形状，落在中间格子的小球数较多，落在两端格子的小球数很少。我们可以增加投球次数，观察小球堆积的分布有无改变；我们也可以改变概率p，观察小球堆积情况的变化；当然，我们也可以增加钉子的数目，看看小球堆积分布的变化情况。模拟Galton钉板试验的步骤如下：（1）确定钉子的位置：将钉子的横、纵坐标存储在两个矩阵X和Y之中。（2）在Galton钉板试验中，小球每次碰到钉子下落时都具有两种可能性。向右的

14、概率为p，向左的概率为q=1-p（整个演示程序中都使用这规则），这里p=0.5，表示向右和向左的机会是相同的。将0,1区间分成两段，区间0,q和p,q+q。如果随机数，让小球向右落下；若，让小球向左落下。将这一过程重复n次，并用直线连接小球落下时所经过的点，这样就模拟了小球从顶端随机地落入某一个格子的过程。（3）模拟小球堆积的形状。输入扔球次数m（例如m50、100、500等），计算落在第格格子的小球数在总球数m中所占的比例，这样当模拟结束时，就得到了频率，用频率反映小球的堆积形状。图5 500个小球堆积频率图main.fig中的主程序G_yanshi.m：function G_yanshi%

15、本函数将完成Calton钉板问题的演示功能load k;load p;m=str2num(k);n=4;y0=2;ballnum=zeros(1,n+1);p=str2num(p);q=1-p;%确定钉子的坐标for i=n+1:-1:1 xding(i,1)=0.5*(n-i);yding(i,1)=(n-i+1)+y0; xliu(i,1)=xding(i,1);yliu(i,1)=yding(i,1)-0.5; for j=2:i xding(i,j)=xding(i,1)+(j-1)*1;yding(i,j)=yding(i,1); xliu(i,j)=xding(i,j);yliu(

16、i,j)=yding(i,j)-0.5; endend%绘制钉板边缘的坐标值x_zuo=2.3000 2.3000 1.8000 1.8000 1.3000 1.3000 0.8000 0.8000 0.3000 0.3000-0.5;y_zuo=6.0000 5.2700 4.7000 4.3000 3.7000 3.3000 2.7000 2.3000 1.7000 1.0000;x_you=2.7000 2.7000 3.2000 3.2000 3.7000 3.7000 4.2000 4.2000 4.7000 4.7000-0.5;y_you=6.0000 5.2700 4.7000

17、 4.3000 3.7000 3.3000 2.7000 2.3000 1.7000 1.0000;%钉子运动过程中的坐标集合xx1=2.5000 2.0000 3.0000 1.5000 2.5000 3.5000 1.0000 2.0000 3.0000 4.0000 0.5000 1.5000 2.5000 3.5000 4.5000-0.5;xx2=2.0000 3.0000 1.5000 2.5000 3.5000 1.0000 2.0000 3.0000 4.0000 0.5000 1.5000 2.5000 3.5000 4.5000-0.5;yy1=5.3000 4.7000

18、4.7000 3.7000 3.7000 3.7000 2.7000 2.7000 2.7000 2.7000 1.7000 1.7000 1.7000 1.7000 1.7000;yy2=4.3000 4.3000 3.3000 3.3000 3.3000 2.3000 2.3000 2.3000 2.3000 1.3000 1.3000 1.3000 1.3000 1.3000;mm=moviein(m); %开始播放动画for h=1:m %绘制边框 for i=2:5 for j=2:5 if xliu(i,j)=0&yliu(i,j)=0 xtt= xliu(i,j), xliu(i

19、,j)+0.3, xliu(i,j)+0.3, xliu(i,j), xliu(i,j)-0.3, xliu(i,j)-0.3,xliu(i,j); ytt=yliu(i,j)+0.5,yliu(i,j)+0.2,yliu(i,j)-0.2,yliu(i,j)-0.5,yliu(i,j)-0.2,yliu(i,j)+0.2,yliu(i,j)+0.5; plot(xtt,ytt) axis(-1,5,0,6); hold on end end end plot(x_zuo,y_zuo,x_you,y_you,b) %以下绘制动态路径 s=rand(1,15); xii=2;yii=6;kk=1

20、;l=1; xi=xx1(1);yi=yy1(1); xs=xx2(1);ys=yy2(1); plot(xii xi,yii yi,g) for j=1:15 if s(j)p l=l+kk; if l15 l=l-kk; end else l=l+kk+1; if l15 l=l-kk-1; end end kk=kk+1; xt=xx1(l);yt=yy1(l); xs=xx2(l-1);ys=yy2(l-1); plot(xi xt xs,yi yt ys, g) xi=xs;yi=ys; end %绘制动态条形图 ballnum(l-10)=ballnum(l-10)+1; ball

21、num1=ballnum./m; bar(0:n,ballnum1) mm(h)=getframe; hold offend save ballnum %保存本次试验的各频数3、动画模拟Galton钉板Poisson分布当二项分布的n很大而p很小时，泊松分布可作为二项分布的近似，其中为np。通常当n10,p0.1时，就可以用泊松公式近似计算。上述公式是泊松分布概率函数，其分布图如图6：图6 泊松分布概率图运行main.fig程序，点击“模拟泊松分布”，根据程序设定，演示结果为“n=10,p=0.1, =1,500个球”既定的结果，如图7：图7 模拟泊松分布概率图从结果图与理论分布图的比较可以看

22、出，演示程序较成功地模拟出了泊松分布。main.fig中泊松分布函数子程序G_bosong.m：function G_bosong%本函数完成泊松演示set(findobj(gcbf,tag,edit0),string,0);set(findobj(gcbf,tag,edit1),string,0);set(findobj(gcbf,tag,edit2),string,0);set(findobj(gcbf,tag,edit3),string,0);set(findobj(gcbf,tag,edit4),string,0);set(findobj(gcbf,tag,edittest),stri

23、ng,0);m=500;n=10;y0=2;ballnum=zeros(1,n+1);p=0.1;q=1-p;%确定钉子的坐标for i=n+1:-1:1 xding(i,1)=0.5*(n-i);yding(i,1)=(n-i+1)+y0; xliu(i,1)=xding(i,1);yliu(i,1)=yding(i,1)-0.5; for j=2:i xding(i,j)=xding(i,1)+(j-1)*1;yding(i,j)=yding(i,1); xliu(i,j)=xding(i,j);yliu(i,j)=yding(i,j)-0.5; endend%绘制钉板边缘的坐标值x_zu

24、o=5.3 5.3 4.8 4.8 4.3 4.3 3.8 3.8 3.3 3.3 2.8 2.8 2.3 2.3 1.8 1.8 1.3 1.3 0.8 0.8 0.3 0.3-0.5;y_zuo=12 11.27 10.7 10.3 9.7 9.3 8.7 8.3 7.7 7.3 6.7 6.3 5.7 5.3 4.7 4.3 3.7 3.3 2.7 2.3 1.7 1;x_you=5.7 5.7 6.2 6.2 6.7 6.7 7.2 7.2 7.7 7.7 8.2 8.2 8.7 8.7 9.2 9.2 9.7 9.7 10.2 10.2 10.7 10.7-0.5;y_you=12

25、11.27 10.7 10.3 9.7 9.3 8.7 8.3 7.7 7.3 6.7 6.3 5.7 5.3 4.7 4.3 3.7 3.3 2.7 2.3 1.7 1;%钉子运动过程中的坐标集合xx1=5.5 5 6 4.5 5.5 6.5 4 5 6 7 3.5 4.5 5.5 6.5 7.5 3 4 5 6 7 8 2.5 3.5 4.5 5.5 6.5 7.5 8.5 2 3 4 5 6 7 8 9 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 1 2 3 4 5 6 7 8 9 10 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5

26、9.5 10.5-0.5;xx2=5 6 4.5 5.5 6.5 4 5 6 7 3.5 4.5 5.5 6.5 7.5 3 4 5 6 7 8 2.5 3.5 4.5 5.5 6.5 7.5 8.5 2 3 4 5 6 7 8 9 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 1 2 3 4 5 6 7 8 9 10 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5-0.5;yy1=linspace(0,0,66);yy1(1)=11.3;for y=2:3; yy1(y)=10.7;endfor y=4:6; yy1(y)=

27、9.7;endfor y=7:10; yy1(y)=8.7;endfor y=11:15; yy1(y)=7.7;endfor y=16:21; yy1(y)=6.7;endfor y=22:28; yy1(y)=5.7;endfor y=29:36; yy1(y)=4.7;endfor y=37:45; yy1(y)=3.7;endfor y=46:55; yy1(y)=2.7;endfor y=56:66; yy1(y)=1.7;endyy2=yy1(2:66)-0.4;mm=moviein(m); %开始播放动画for h=1:m %绘制边框 for i=2:11 for j=2:11

28、if xliu(i,j)=0&yliu(i,j)=0 xtt= xliu(i,j), xliu(i,j)+0.3, xliu(i,j)+0.3, xliu(i,j), xliu(i,j)-0.3, xliu(i,j)-0.3,xliu(i,j); ytt=yliu(i,j)+0.5,yliu(i,j)+0.2,yliu(i,j)-0.2,yliu(i,j)-0.5,yliu(i,j)-0.2,yliu(i,j)+0.2,yliu(i,j)+0.5; plot(xtt,ytt) axis(-1,11,0,12); hold on end end end plot(x_zuo,y_zuo,x_yo

29、u,y_you,b) %以下绘制动态路径 s=rand(1,66); xii=5;yii=12;kk=1;l=1; xi=xx1(1);yi=yy1(1); xs=xx2(1);ys=yy2(1); plot(xii xi,yii yi,g) for j=1:66 if s(j)p l=l+kk; if l66 l=l-kk; end else l=l+kk+1; if l66 l=l-kk-1; end end kk=kk+1; xt=xx1(l);yt=yy1(l); xs=xx2(l-1);ys=yy2(l-1); plot(xi xt xs,yi yt ys, g) xi=xs;yi=

30、ys; end %绘制动态条形图 ballnum(l-55)=ballnum(l-55)+1; ballnum1=ballnum./m; bar(0:n,ballnum1) mm(h)=getframe; hold offend save ballnum %保存本次试验的各频数G_sum4、动画模拟Galton钉板几何分布几何分布（Geometric distribution）是离散型概率分布。其中一种定义为：在第n次伯努利试验中，试验k次才得到第一次成功的机率。详细的说，是：前k-1次皆失败，第k次成功的概率。其公式为：我们在演示程序中，用500个球模拟k=10,p=0.1时，小球落入1号格

31、子的频数。通过计算，我们得到理论频数为：0.99*0.1*50019.3，而我们演示结果如图8：图8 模拟几何分布结果有一点要解释的是，因为我们要统计的是小球落入1号格子的频数，所以对其他落入其他格子的小球数没有进行统计。从结果得出模拟频数为17，与理论值19.3差距不大，说明我们的演示能较成功地模拟出几何分布。main.fig中几何分布函数子程序G_jihe.m：function G_jihe%本函数完成几何分布演示set(findobj(gcbf,tag,edit0),string,0);set(findobj(gcbf,tag,edit1),string,0);set(findobj(g

32、cbf,tag,edit2),string,0);set(findobj(gcbf,tag,edit3),string,0);set(findobj(gcbf,tag,edit4),string,0);set(findobj(gcbf,tag,edittest),string,0);m=500;n=10;y0=2;ballnum=zeros(1,n+1);p=0.1;q=1-p;pd2=0;pd=0;%确定钉子的坐标for i=n+1:-1:1 xding(i,1)=0.5*(n-i);yding(i,1)=(n-i+1)+y0; xliu(i,1)=xding(i,1);yliu(i,1)=

33、yding(i,1)-0.5; for j=2:i xding(i,j)=xding(i,1)+(j-1)*1;yding(i,j)=yding(i,1); xliu(i,j)=xding(i,j);yliu(i,j)=yding(i,j)-0.5; endend%绘制钉板边缘的坐标值x_zuo=5.3 5.3 4.8 4.8 4.3 4.3 3.8 3.8 3.3 3.3 2.8 2.8 2.3 2.3 1.8 1.8 1.3 1.3 0.8 0.8 0.3 0.3-0.5;y_zuo=12 11.27 10.7 10.3 9.7 9.3 8.7 8.3 7.7 7.3 6.7 6.3 5.

34、7 5.3 4.7 4.3 3.7 3.3 2.7 2.3 1.7 1;x_you=5.7 5.7 6.2 6.2 6.7 6.7 7.2 7.2 7.7 7.7 8.2 8.2 8.7 8.7 9.2 9.2 9.7 9.7 10.2 10.2 10.7 10.7-0.5;y_you=12 11.27 10.7 10.3 9.7 9.3 8.7 8.3 7.7 7.3 6.7 6.3 5.7 5.3 4.7 4.3 3.7 3.3 2.7 2.3 1.7 1;%钉子运动过程中的坐标集合xx1=5.5 5 6 4.5 5.5 6.5 4 5 6 7 3.5 4.5 5.5 6.5 7.5 3

35、4 5 6 7 8 2.5 3.5 4.5 5.5 6.5 7.5 8.5 2 3 4 5 6 7 8 9 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 1 2 3 4 5 6 7 8 9 10 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5-0.5;xx2=5 6 4.5 5.5 6.5 4 5 6 7 3.5 4.5 5.5 6.5 7.5 3 4 5 6 7 8 2.5 3.5 4.5 5.5 6.5 7.5 8.5 2 3 4 5 6 7 8 9 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5

36、1 2 3 4 5 6 7 8 9 10 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5-0.5;yy1=linspace(0,0,66);yy1(1)=11.3;for y=2:3; yy1(y)=10.7;endfor y=4:6; yy1(y)=9.7;endfor y=7:10; yy1(y)=8.7;endfor y=11:15; yy1(y)=7.7;endfor y=16:21; yy1(y)=6.7;endfor y=22:28; yy1(y)=5.7;endfor y=29:36; yy1(y)=4.7;endfor y=37:45;

37、 yy1(y)=3.7;endfor y=46:55; yy1(y)=2.7;endfor y=56:66; yy1(y)=1.7;endyy2=yy1(2:66)-0.4;mm=moviein(m); %开始播放动画for h=1:m %绘制边框 for i=2:11 for j=2:11 if xliu(i,j)=0&yliu(i,j)=0 xtt= xliu(i,j), xliu(i,j)+0.3, xliu(i,j)+0.3, xliu(i,j), xliu(i,j)-0.3, xliu(i,j)-0.3,xliu(i,j); ytt=yliu(i,j)+0.5,yliu(i,j)+0.2,yliu(i,j)-0.2,yliu(i,j)-0.5,yliu(i,j)-0.2,yliu(i,j)+0.2,yliu(i,j)