哈工大数学实验大作业 Galton钉板问题Word文档下载推荐.docx

1、问题1、小球自上方落下，经过n个钉子。每经过一个钉子时只有两种可能结果：向右或向左。这是一个具有两个结果（成功和失败）的随机试验E，将向右视为成功，成功的概率为p，向左为失败，失败的概率为q=1-p。小球碰到一个钉子下落一格，相当于进行了一次试验E。小球自顶端落下，碰到n个钉子，最终落在某个格子的过程，恰好相当于将试验E重复了n次，因此一次投球过程就是一个n重贝努利试验（将仅有两个相互排斥结果的试验E独立重复n次，构成了n重贝努利试验En）。n重贝努利试验的成功次数X正好是小球向右移动的次数，它是一个随机变量。根据概率论的结果有XB（n,p）。对于一个随机变量，我们首先要弄清楚它的取值范围，X

2、的取值范围为0,1,2,n，这是什么意思呢？在Galton钉板模型中X0表示小球向右移动的次数，也就是小球一直向左移动，所以它恰好要落在编号为0的格子里；同理X1表示小球恰好要落在编号为1的格子里，依次类推，这就是说，X是小球最终落进的格子编号数，当然它也对应为小球向右移动的次数。二项随机变量的分布列为： 则在Galton钉板实验中，p0.5，底槽中各格的频数应为k*pi。问题2、演示程序已经完成，具体程序会在下文“实验过程”中体现，以下是演示结果如图1：图1 Galton钉板问题演示结果可以从界面看到，我们的演示程序能方便地对k和p的值进行修改再演示，同时还能统计落入每个格子的小球数。在演示

3、界面下方，有一列是用来检验2上侧分位数是否满足等式P22（n）=我们在程序中设置了=（0.005,0.01,0.025,0.05,0.10,0.25,0.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模拟几何分布演示结果鉴于模拟的难度，我们在演示程序中设定了参数不可更改，尽管可能导致结果有偏差，但我们经

4、过多次演示，发现偏差极小。四、背景知识介绍1、随机变量随机变量是随机试验结果的函数，其特点是在试验前，并不能预知这个函数将取何值，这要凭机会，就是“随机”的意思。一旦试验后，取值就确定了。例如，我在3月31日买了一张奖券，到6月30日开奖。当我买这张奖券时，有人可以对我说：“你中奖的金额是个随机的变量，其值在6月30日抽奖试验做过之后才能确定。”明白了这一点就不难举出许多随机变量的例子。例如，某出租车公司的电话订车中心，一天内接到的订车电话的次数；某射手对一活动靶进行射击，到击中目标为止，所进行的射击次数；从一批灯泡中，任取一只，测定这只灯泡的寿命，等等，这些都是随机变量。2、n重Bernou

5、lli试验当依照一定的质量标准，从大批产品中抽出一件进行产品质量合格性检查时，得出的结果可以是二者之一：“这是件不合格产品”或“这是件合格产品”。如果将这样的抽检一件产品看作是进行一次试验，则试验的结果可以是发生A（这是件不合格产品）或（A不发生，即产品质量合格）。称这种只有两个可能结果A（称“成功”）或（称“失败”）的试验为Bernoulli试验。有很多试验，其可能的结果不止两个，但由于人们常只对试验是否发生某一种特定结果感兴趣，因而可将之归结为Bernoulli试验、例如，明天的天气可以有多种情况，但若只关心明天是否下雨，则观察明天的天气（作为一次试验），其结果只有两个：“下雨”或“不下雨

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

7、均每小时有3分钟要使用外线，在考虑该总机应设置多少条外线合适的问题时，可归结为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试验模型，这里n5

8、000。3、二项分布在“成功”概率是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）如果想保持结果一致，可与ran

9、d（seed,s）配合使用，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钉板试验运行观察程序mai

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

11、的概率为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

12、%本函数将完成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

13、（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.700

14、0 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

15、 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 if xliu（i,j）=0&yliu（i,j）=0 xtt= xliu（i,j）, xliu（i,j）+0.3,

16、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 plot（x_zuo,y_zuo,x_you,y_you,b） %以下绘制动态路径 s=rand（1,15）; xii=2;yii=6;kk=1;l=1; xi=xx1（1）;yi=yy

17、1（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 l l=l-kk; else l=l+kk+1; l=l-kk-1; kk=kk+1; xt=xx1（l）;yt=yy1（l）; xs=xx2（l-1）;ys=yy2（l-1）; plot（xi xt xs,yi yt ys, xi=xs;yi=ys; %绘制动态条形图 ballnum（l-10）=ballnum（l-10）+1; ballnum1=ballnum./m; bar（0:n,ballnum1） mm（h）=getframe;

18、 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 模拟泊松分布概率图从结果图与理论分布图的比较可以看出，演示程序较成功地模拟出了泊松分布。main.fig中泊松分布函数子程序G_bosong.m：fu

19、nction G_bosong%本函数完成泊松演示set（findobj（gcbf,tag,edit0）,string0）;edit1edit2edit3edit4edittestm=500;n=10;p=0.1;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.7 5.3 4.7 4.3 3.7 3.3 2.7 2.3 1.7 1;

20、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 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

21、 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 1 2 3 4 5 6 7 8 9 10 0.5 1.5 2.5 3.5 4.5 5.5 6.5

22、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;for y=4:6; yy1（y）=9.7;for y=7:10; yy1（y）=8.7;for y=11:15; yy1（y）=7.7;for y=16:21; yy1（y）=6.7;for y=22:28; yy1（y）=5.7;for y=29:36; yy1（y）=4.7;for y=37:45; yy1（y）=3.7;for y=46:55; yy1（y）=2.7;for y=56:66; yy1（y）=1.7;yy2=yy1（2:

23、66）-0.4;11 axis（-1,11,0,12）; s=rand（1,66）; xii=5;yii=12;66 ballnum（l-55）=ballnum（l-55）+1;G_sum4、动画模拟Galton钉板几何分布几何分布（Geometric distribution）是离散型概率分布。其中一种定义为：在第n次伯努利试验中，试验k次才得到第一次成功的机率。详细的说，是：前k-1次皆失败，第k次成功的概率。其公式为：我们在演示程序中，用500个球模拟k=10,p=0.1时，小球落入1号格子的频数。通过计算，我们得到理论频数为：0.99*0.1*50019.3，而我们演示结果如图8：图8 模拟几何分布结果有一点要解释的是，因为我们要统计的是小球落入1号格子的频数，所以对其他落入其他格子的小球数没有进行统计。从结果得出模拟频数为17，与理论值19.3差距不大，说明我们的演示能较成功地模拟出几何分布。main.fig中几何分布函数子程序G_jihe.m：function G_jihe%本函数完成几何分布演示set（