产品
摘要
在《概率论与数理统计》学习中，阐释过中心极限定理，并对中心极限定理的应用形式有过介绍。在其中，高尔顿钉板实验常常作为应用实例，利用中心极限定理来进行解释。但是课上展示环节主要利用几何画板或者利用实际手工操作实验视频来进行展示，有局限性，不能模拟多种情况，比如钉子层数的不同、实验小球的个数、实验开始时第一层钉子个数对于实验结果的影响等等。本文主要基于Matlab软件，对高尔顿实验进行模拟，从而实现较为直观的对高尔顿实验的多种情况进行模拟。
关键词：中心极限定理 高尔顿钉板实验 Matlab
1、 问题描述
高尔顿钉板实验，常作为教学实验在《概率论与数理统计》中进行展示。其设计者为英国生物统计学家高尔顿。
1.1 实验装置
由钉子、略小于钉子间隔的小球、带有各自的底板等组成，如图表 1所示，图表来源见水印。
图表 1高尔顿实验装置模型
1.2 实验过程
从入口处放进一个直径略小于两颗钉子之间的距离的小圆玻璃球，小圆球由于重力向下运动，运动过程中碰到钉子，碰到钉子后向右或者向左的概率相等，都为1/2。以此类推，直到滚到底板的一个格子内为止。把实验设定数量且同样大小的小球不断从入口处放下，只要球的数目相当大，它们在底板将堆成近似于正态分布的密度函数图形（即：中间高，两边低，左右对称），其中n为钉子的层数。实验操作如视频 1 高尔顿实验操作视频
2、 基本原理
2.1 中心极限定理
2.1.1 林德贝尔格-勒维定理
具有有限方差的一列同分布的随机变量的和经过标准化后是以标准正态分布为极限的，这就是独立同分布的中心极限定理或称为林德贝尔格-勒维中心极限定理。
设1）{Xn} 独立同分布；2）∀k, E(Xk)= , D(Xk) = 2,0 < 2 <∞，这时，E()=n, D()=n2
则∀x∈R，=，等价于=
应用形式（n较大时）：~N(0,1)或~N()
~N(0,1)或~N()
2.1.2 棣莫弗-拉普拉斯中心极限定理
设 n=1,2,3…(0<p<1),这时E=np, D=np(1-p)
则∀x∈R，=
应用形式：X~B(n,p),n较大,0<p<1
有~N(0,1)或X~N(np,np(1-p))
2.2 高尔顿实验解释
实验中挡板钉子层数为n，下落的小球每碰到一个钉子后向左走的概率为p，则向右走的概率为1-p；小球从顶端槽进入装置，每个小球中间总共向左走的步数为n1，则向右走的步数为n-n1记为n2，落入槽的位置为N，不妨令没有钉子小球垂直落到的位置为N=n，则有：
N=n+(n1-n2)=2n1→n1=N/2(n1=0,1,2,3…)
因此，实验中我们确定槽的个数为n,而把槽的位置从左至右依次编号为0，2，4，6,…,2n，这样做的目的是使实验便于分析。根据中心极限定理采用二项式概率分布，一个小球落下，出现左走步数n1，右走步数n2的概率为：
X~N(p,1-p)
P(n1)=P(N)=P(N/2)=
这也是小球落入位置的概率分布，是一个正态二项分布。
3、 Matlab模拟
m=100,n=6,y0=3,w=10000,v=1000;设置需要模拟实验的钉子层数n、作图开始纵坐标y0、小球个数等
ballnum=zeros(1,n+1); 生成1行n+1列零矩阵
p=0.5,q=1-p; 设置小球向左向右下落的概率
for i=n+1:-1:1 创建钉子的坐标x,y
x(i,1)=0.5*(n-i+1); 得到钉子的横坐标，为直观显示正太分布结果（即将之后条形图的范围设置在1,2,3…附近，因此将钉子的很坐标设为0.5,1.5,2.5,3.5,…）
y(i,1)=(n-i+1)+y0; 得到钉子的纵坐标，因要在钉子坐标下画出最后的条形图， 故不能y坐标的起始值设为0，设为y0。
for j=2:1:i
x(i,j)=x(i,1)+(j-1)*1;
y(i,j)=y(i,1);
end
end
mm=moviein(m) 调用moviein函数对内存进行初始化，创建一个足够大的矩阵，使之能够容纳基于当前坐标轴大小的一系列指定的图形（此处称为帧）。
for i=1:1:m
s=rand(1,w); 产生w个随机数
xi=x(1,1); 小球遇到第一个钉子
yi=y(1,1);
k=1;d=1;
for j=1:1:n
plot(x(1:n,:),y(1:n,:),'o',x(n+1,:),y(n+1,:),'-')画钉子的位置
axis([-2 n+2 0 y0+n+1]),hold on
k=k+1; 小球下落一格
if s(j)>p
d=d+0; 小球左移
else
d=d+1; 小球右移
end
xt=x(k,d);yt=y(k,d); 小球下落点的坐标
h=plot([xi,xt],[yi,yt]);axis([-2 n+2 0 y0+n+1])；画小球的运动轨迹
xi=xt;yi=yt;
end
ballnum(d)=ballnum(d)+1; 计数
ballnum1=3*ballnum./m;
title('高尔顿matlab实验模拟');
bar([0:n],ballnum1),axis([-2 n+2 0 y0+n+1])画各格子的频率
mm(i)=getframe; 调用getframe函数生成每个帧。该函数返回一个列矢量，利用这个矢量，就可以创建一个电影动画矩阵。
hold off
end
完整matlab程序：
m=100,n=6,y0=3,w=10000,v=1000
ballnum=zeros(1,n+1);
p=0.5,q=1-p
for i=n+1:-1:1
x(i,1)=0.5*(n-i+1);
y(i,1)=(n-i+1)+y0;
for j=2:1:i
x(i,j)=x(i,1)+(j-1)*1;
y(i,j)=y(i,1);
end
end
mm=moviein(m)
for i=1:1:m
s=rand(1,w);
xi=x(1,1);
yi=y(1,1);
k=1;d=1;
for j=1:1:n
plot(x(1:n,:),y(1:n,:),'o',x(n+1,:),y(n+1,:),'-')
axis([-2 n+2 0 y0+n+1]),hold on
k=k+1;
if s(j)>p
d=d+0;
else
d=d+1;
end
xt=x(k,d);yt=y(k,d);
h=plot([xi,xt],[yi,yt]);axis([-2 n+2 0 y0+n+1])
xi=xt;yi=yt;
end
ballnum(d)=ballnum(d)+1;
ballnum1=3*ballnum./m;
title('高尔顿matlab实验模拟');
bar([0:n],ballnum1),axis([-2 n+2 0 y0+n+1])
mm(i)=getframe;
hold off
end
4、 模拟结果
如视频 2所示： 高尔顿钉板实验matlab模拟
5、 结论与心得体会
利用zeros函数生成零矩阵；利用rand函数随机生成矩阵；利用for循环语句，生成钉子的坐标；利用plot函数画出钉子坐标以及利用bar函数生成条形图画出最后的模拟结果。利用moviein以及getframe函数捕捉每一个帧的画面并将数据保存，在循环结束后进行播放，更加直观生动。建立function函数，只需要将模拟层数以及钉子数目等输入其中，就可以模拟多种情况。
高尔顿实验在n较大时，更为接近正态分布曲线，同时也利用高尔顿实验验证了中心极限定理。
参考文献
[1]侯臣平,矫媛媛.Matlab在《概率论与数理统计》教学中的应用[J].教育教学论坛,2019(05):156-157. [2]武新乾,张翠霞,轩凤霞.高尔顿钉板试验动态图形软件的设计与制作[J].洛阳师范学院学报,2015,34(02):75-78.
155-2731-8020
