哈工大数学实验大作业 Galton钉板问题.docx
- 文档编号:3924779
- 上传时间:2022-11-26
- 格式:DOCX
- 页数:20
- 大小:173.54KB
哈工大数学实验大作业 Galton钉板问题.docx
《哈工大数学实验大作业 Galton钉板问题.docx》由会员分享,可在线阅读,更多相关《哈工大数学实验大作业 Galton钉板问题.docx(20页珍藏版)》请在冰豆网上搜索。
哈工大数学实验大作业Galton钉板问题
Galton钉板问题
一、实验目的与要求
1.复习概率论中随机变量、概率分布、二项分布、均值和分布函数等概念。
2.理解Galton钉板实验中小球落入格子所服从的规律。
3.了解Matlab软件中进行动画演示的命令。
4.掌握Matlab软件中进行随机模拟的方法。
二、问题描述
所有现象的“因”和“果”,即“条件”和“结果”之间在客观上都存在着一定的规律,这种规律通常可以分成两类:
一是确定性的规律,另一类是非确定性规律。
对于确定性的系统,当已知条件是充分时,那么实验的结果也是确定的,即在每一次试验进行以前,可以预见试验产生的结果。
但若条件不充分时,就无法预测试验的结果,这就产生了“因果律的缺失”的随机现象。
随机现象在实践中是大量遇到的,如掷骰子。
虽然无法由“因”预测“果”,但是当进行大量重复试验时,因果之间仍会呈现一种统计规律。
概率方法建立在“重复试验”的基础上,统计规律只有在大量重复后才会呈现出来,诸如随机变量。
分布、均值、方差等概念无一不体现了重复的思想。
以下围绕着Galton钉板模型来讨论。
Galton钉板试验是由英国生物统计学家Galton设计的。
在一板上有n排钉子,图0所示的是n=4的情况。
图中10个圆点表示10颗钉子,在钉子的下方有5个格子,分别编号为0,1,2,3,4。
自Galton钉板的上方扔进一个小球任其自由下落,在下落的过程中当小球碰到钉子时,从左边落下与从右边落下的机会相等、碰到下一排钉子时又是如此。
最后落入底板中的某一个格子。
图0Galton钉板模型(n=4)
问题:
1、一个一个从顶部放入k个小球,底槽中各格的理论频数应为多少?
2、
(1)编写一演示程序(用MATLAB语言),其界面与完成的基本功能如书中P121.图8.9所示;
(2)利用自由度为4的Χ2分布临界值表及统计量
(其中vi为0~4号格子中第i个实验频数)。
检验假设
H0:
所实验的分布服从B(4,p)
3、用类似方法模拟Poisson分布与几何分布。
三、问题分析
问题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
forj=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
end
end
end
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>15
l=l-kk;
end
else
l=l+kk+1;
ifl>15
l=l-kk-1;
end
end
kk=kk+1;
xt=xx1(l);yt=yy1(l);
xs=xx2(l-1);ys=yy2(l-1);
plot([xixtxs],[yiytys],'g')
xi=xs;yi=ys;
end
%绘制动态条形图
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');
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'),'string','0');
m=500;n=10;y0=2;
ballnum=zeros(1,n+1);
p=0.1;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=[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;end
fory=4:
6;yy1(y)=9.7;end
fory=7:
10;yy1(y)=8.7;end
fory=11:
15;yy1(y)=7.7;end
fory=16:
21;yy1(y)=6.7;end
fory=22:
28;yy1(y)=5.7;end
fory=29:
36;yy1(y)=4.7;end
fory=37:
45;yy1(y)=3.7;end
fory=46:
55;yy1(y)=2.7;end
fory=56:
66;yy1(y)=1.7;end
yy2=yy1(2:
66)-0.4;
mm=moviein(m);%开始播放动画
forh=1:
m
%绘制边框
fori=2:
11
forj=2:
11
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,11,0,12]);
holdon
end
end
end
plot(x_zuo,y_zuo,x_you,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([xiixi],[yiiyi],'g')
forj=1:
66
ifs(j)>p
l=l+kk;
ifl>66
l=l-kk;
end
else
l=l+kk+1;
ifl>66
l=l-kk-1;
end
end
kk=kk+1;
xt=xx1(l);yt=yy1(l);
xs=xx2(l-1);ys=yy2(l-1);
plot([xixtxs],[yiytys],'g')
xi=xs;yi=ys;
end
%绘制动态条形图
ballnum(l-55)=ballnum(l-55)+1;
ballnum1=ballnum./m;
bar([0:
n],ballnum1)
mm(h)=getframe;
holdoff
end
saveballnum%保存本次试验的各频数
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(findobj(gcbf,'tag','edit0'),'string','0');
set(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 哈工大数学实验大作业 Galton钉板问题 哈工大 数学 实验 作业 Galton 问题