哈工大数学实验大作业 Galton钉板问题Word文档下载推荐.docx
《哈工大数学实验大作业 Galton钉板问题Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《哈工大数学实验大作业 Galton钉板问题Word文档下载推荐.docx（20页珍藏版）》请在冰豆网上搜索。
问题1、小球自上方落下，经过n个钉子。
每经过一个钉子时只有两种可能结果：
向右或向左。
这是一个具有两个结果（成功和失败）的随机试验E，将向右视为成功，成功的概率为p，向左为失败，失败的概率为q=1-p。
小球碰到一个钉子下落一格，相当于进行了一次试验E。
小球自顶端落下，碰到n个钉子，最终落在某个格子的过程，恰好相当于将试验E重复了n次，因此一次投球过程就是一个n重贝努利试验（将仅有两个相互排斥结果的试验E独立重复n次，构成了n重贝努利试验En）。
n重贝努利试验的成功次数X正好是小球向右移动的次数，它是一个随机变量。
根据概率论的结果有X~B（n,p）。
对于一个随机变量，我们首先要弄清楚它的取值范围，X的取值范围为0,1,2…,n，这是什么意思呢？
在Galton钉板模型中X＝0表示小球向右移动的次数，也就是小球一直向左移动，所以它恰好要落在编号为0的格子里；
同理X＝1表示小球恰好要落在编号为1的格子里，依次类推，这就是说，X是小球最终落进的格子编号数，当然它也对应为小球向右移动的次数。
二项随机变量
的分布列为：
则在Galton钉板实验中，p＝0.5，底槽中各格的频数应为k*pi。
问题2、演示程序已经完成，具体程序会在下文“实验过程”中体现，以下是演示结果如图1：
图1Galton钉板问题演示结果
可以从界面看到，我们的演示程序能方便地对k和p的值进行修改再演示，同时还能统计落入每个格子的小球数。
在演示界面下方，有一列是用来检验Χ2上侧α分位数是否满足等式
P[Χ2≥Χα2（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
图2HELP文档
问题3、Poisson分布与几何分布演示结果如图3和图4：
图3Galton模拟Poisson分布演示结果
图4Galton模拟几何分布演示结果
鉴于模拟的难度，我们在演示程序中设定了参数不可更改，尽管可能导致结果有偏差，但我们经过多次演示，发现偏差极小。
四、背景知识介绍
1、随机变量
随机变量是随机试验结果的函数，其特点是在试验前，并不能预知这个函数将取何值，这要凭机会，就是“随机”的意思。
一旦试验后，取值就确定了。
例如，我在3月31日买了一张奖券，到6月30日开奖。
当我买这张奖券时，有人可以对我说：
“你中奖的金额ξ是个随机的变量，其值在6月30日‘抽奖试验’做过之后才能确定。
”
明白了这一点就不难举出许多随机变量的例子。
例如，某出租车公司的电话订车中心，一天内接到的订车电话的次数ξ；
某射手对一活动靶进行射击，到击中目标为止，所进行的射击次数η；
从一批灯泡中，任取一只，测定这只灯泡的寿命ζ，等等，这些都是随机变量。
2、n重Bernoulli试验
当依照一定的质量标准，从大批产品中抽出一件进行产品质量合格性检查时，得出的结果可以是二者之一：
“这是件不合格产品”或“这是件合格产品”。
如果将这样的抽检一件产品看作是进行一次试验，则试验的结果可以是发生A（这是件不合格产品）或Ā（A不发生，即产品质量合格）。
称这种只有两个可能结果A（称“成功”）或Ā（称“失败”）的试验为Bernoulli试验。
有很多试验，其可能的结果不止两个，但由于人们常只对试验是否发生某一种特定结果感兴趣，因而可将之归结为Bernoulli试验、例如，明天的天气可以有多种情况，但若只关心明天是否下雨，则观察明天的天气（作为一次试验），其结果只有两个：
“下雨”或“不下雨”，因而可被看作是个Bernoulli试验。
实际上常要考察独立重复进行Bernoulli试验的序列，并将这一独立重复试验序列作为单独一个复合试验来对待。
这里，所谓独立重复进行Bernoulli试验的意思是，这个序列中的每一试验的结果都只能发生A或Ā，且发生A的概率一样，是某个值p=P（A）。
当然，发生Ā的概率也就是一样地是q=P（Ā），并且每一试验发生的结果不会影响其他试验出现的结果。
独立重复试验序列最重要的特性是序列由独立重复进行n次Bernoulli试验组成，简称为n重Bernoulli试验。
在很多问题中可以用上n重Bernoulli试验模型。
例如，若学校的电话总机设有99个分机，已知每号分机平均每小时有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试验模型，这里n＝5000。
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）
产生m×
n个（0,1）区间中的随机数，并将这些随机数存于一个m×
n矩阵中。
每次调用rand（m,n）的结果都会不同。
rand（‘seed’,s）
如果想保持结果一致，可与rand（‘seed’,s）配合使用，s为一正整数。
Moviein（n）
创建动画矩阵，制作动画矩阵数据
Getframe
拷贝动画矩阵
Movie（mat,m）
播放动画矩阵m次
例如：
输入命令：
rand（‘seed’,1）,u=rand（1,6）
得到结果：
0.51290.46050.35040.09500.43370.7092
而且再次运行该指令时结果保持不变，除非重新设置种子seed的值。
rand（‘seed’,2）,u=rand（1,6）
•得到结果：
u=0.02580.92100.70080.19010.86730.4185
这样结果才会产生变化。
2、动画模拟Galton钉板试验
运行观察程序main.fig，屏幕将出现图形窗口，动画模拟扔球过程。
如图5是模拟向一个4层Galton钉板扔500次小球的过程的最后的结果。
在模拟过程中我们看到，每一个小球落在哪一个格子是无法预测的，但小球逐渐堆积成一种单峰的形状，落在中间格子的小球数较多，落在两端格子的小球数很少。
我们可以增加投球次数，观察小球堆积的分布有无改变；
我们也可以改变概率p，观察小球堆积情况的变化；
当然，我们也可以增加钉子的数目，看看小球堆积分布的变化情况。
模拟Galton钉板试验的步骤如下：
（1）确定钉子的位置：
将钉子的横、纵坐标存储在两个矩阵X和Y之中。
（2）在Galton钉板试验中，小球每次碰到钉子下落时都具有两种可能性。
向右的概率为p，向左的概率为q=1-p（整个演示程序中都使用这规则），这里p=0.5，表示向右和向左的机会是相同的。
将[0,1]区间分成两段，区间[0,q]和[p,q+q]。
如果随机数
，让小球向右落下；
若
，让小球向左落下。
将这一过程重复n次，并用直线连接小球落下时所经过的点，这样就模拟了小球从顶端随机地落入某一个格子的过程。
（3）模拟小球堆积的形状。
输入扔球次数m（例如m＝50、100、500等），计算落在第
格格子的小球数在总球数m中所占的比例，这样当模拟结束时，就得到了频率
，用频率反映小球的堆积形状。
图5500个小球堆积频率图
main.fig中的主程序G_yanshi.m：
functionG_yanshi
%本函数将完成Calton钉板问题的演示功能
loadk;
loadp;
m=str2num（k）;
n=4;
y0=2;
ballnum=zeros（1,n+1）;
p=str2num（p）;
q=1-p;
%确定钉子的坐标
fori=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;
forj=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;
end
end
%绘制钉板边缘的坐标值
x_zuo=[2.30002.30001.80001.80001.30001.30000.80000.80000.30000.3000]-0.5;
y_zuo=[6.00005.27004.70004.30003.70003.30002.70002.30001.70001.0000];
x_you=[2.70002.70003.20003.20003.70003.70004.20004.20004.70004.7000]-0.5;
y_you=[6.00005.27004.70004.30003.70003.30002.70002.30001.70001.0000];
%钉子运动过程中的坐标集合
xx1=[2.50002.00003.00001.50002.50003.50001.00002.00003.00004.00000.50001.50002.50003.50004.5000]-0.5;
xx2=[2.00003.00001.50002.50003.50001.00002.00003.00004.00000.50001.50002.50003.50004.5000]-0.5;
yy1=[5.30004.70004.70003.70003.70003.70002.70002.70002.70002.70001.70001.70001.70001.70001.7000];
yy2=[4.30004.30003.30003.30003.30002.30002.30002.30002.30001.30001.30001.30001.30001.3000];
mm=moviein（m）;
%开始播放动画
forh=1:
m
%绘制边框
fori=2:
5
ifxliu（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,5,0,6]）;
holdon
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=yy1
（1）;
xs=xx2
（1）;
ys=yy2
（1）;
plot（[xiixi],[yiiyi],'
g'
forj=1:
15
ifs（j）>
p
l=l+kk;
ifl>
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（[xixtxs],[yiytys],'
xi=xs;
yi=ys;
%绘制动态条形图
ballnum（l-10）=ballnum（l-10）+1;
ballnum1=ballnum./m;
bar（[0:
n],ballnum1）
mm（h）=getframe;
holdoff
end
saveballnum%保存本次试验的各频数
3、动画模拟Galton钉板Poisson分布
当二项分布的n很大而p很小时，泊松分布可作为二项分布的近似，其中λ为np。
通常当n≧10,p≦0.1时，就可以用泊松公式近似计算。
上述公式是泊松分布概率函数，其分布图如图6：
图6泊松分布概率图
运行main.fig程序，点击“模拟泊松分布”，根据程序设定，演示结果为“n=10,p=0.1,λ=1,500个球”既定的结果，如图7：
图7模拟泊松分布概率图
从结果图与理论分布图的比较可以看出，演示程序较成功地模拟出了泊松分布。
main.fig中泊松分布函数子程序G_bosong.m：
functionG_bosong
%本函数完成泊松演示
set（findobj（gcbf,'
tag'
'
edit0'
）,'
string'
0'
）;
edit1'
edit2'
edit3'
edit4'
edittest'
m=500;
n=10;
p=0.1;
x_zuo=[5.35.34.84.84.34.33.83.83.33.32.82.82.32.31.81.81.31.30.80.80.30.3]-0.5;
y_zuo=[1211.2710.710.39.79.38.78.37.77.36.76.35.75.34.74.33.73.32.72.31.71];
x_you=[5.75.76.26.26.76.77.27.27.77.78.28.28.78.79.29.29.79.710.210.210.710.7]-0.5;
y_you=[1211.2710.710.39.79.38.78.37.77.36.76.35.75.34.74.33.73.32.72.31.71];
xx1=[5.5564.55.56.545673.54.55.56.57.53456782.53.54.55.56.57.58.5234567891.52.53.54.55.56.57.58.59.5123456789100.51.52.53.54.55.56.57.58.59.510.5]-0.5;
xx2=[564.55.56.545673.54.55.56.57.53456782.53.54.55.56.57.58.5234567891.52.53.54.55.56.57.58.59.5123456789100.51.52.53.54.55.56.57.58.59.510.5]-0.5;
yy1=linspace（0,0,66）;
yy1
（1）=11.3;
fory=2:
3;
yy1（y）=10.7;
fory=4:
6;
yy1（y）=9.7;
fory=7:
10;
yy1（y）=8.7;
fory=11:
15;
yy1（y）=7.7;
fory=16:
21;
yy1（y）=6.7;
fory=22:
28;
yy1（y）=5.7;
fory=29:
36;
yy1（y）=4.7;
fory=37:
45;
yy1（y）=3.7;
fory=46:
55;
yy1（y）=2.7;
fory=56:
66;
yy1（y）=1.7;
yy2=yy1（2:
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_sum
4、动画模拟Galton钉板几何分布
几何分布（Geometricdistribution）是离散型概率分布。
其中一种定义为：
在第n次伯努利试验中，试验k次才得到第一次成功的机率。
详细的说，是：
前k-1次皆失败，第k次成功的概率。
其公式为：
我们在演示程序中，用500个球模拟k=10,p=0.1时，小球落入1号格子的频数。
通过计算，我们得到理论频数为：
0.9^9*0.1*500≈19.3，而我们演示结果如图8：
图8模拟几何分布结果
有一点要解释的是，因为我们要统计的是小球落入1号格子的频数，所以对其他落入其他格子的小球数没有进行统计。
从结果得出模拟频数为17，与理论值19.3差距不大，说明我们的演示能较成功地模拟出几何分布。
main.fig中几何分布函数子程序G_jihe.m：
functionG_jihe
%本函数完成几何分布演示
set（
升级会员 