1701009700
①确定钉子的位置。将钉子的横、纵坐标存储在一个矩阵中；
1701009701
1701009702
1701009703
②模拟了小球从顶端随机地落入某一格子的过程。设向右的概率为p，向左的概率为q=1-p；将[0,1]分成两段，区间[p,p]和(p,1]。利用rand[ ]产生一个介于0和1之间的随机数u，如果随机数，让小球落向左边，否则落向右边；将这一过程重复n次，并用直线连接小球落下时所经过的点。
1701009704
1701009705
1701009706
③模拟小球堆积的形状。输入扔球次数m，计算落在第i个格子的小球数mi在总球数m中所占的比例，这样当模拟结束时，就得到了频率，用频率反映小球堆积的形状。
1701009707
1701009708
④利用movie完成实时动画设计，观察每一次实验结果。
1701009709
1701009710
选取p=0.5，编写MATLAB程序如下：
1701009711
1701009712
clc,clear,close all %清屏和清除变量 warning off %消除警告 %Galton钉板试验 %设置参数 m = 100; n=6; y0=2; ballnum = zeros(1,n+1); %初始化设置 p=0.5; %概率 q=1-p; figure(‘color’,[1,1,1]) %设置图形背景为白色 for i=n+1:-1:1 %创建钉子的坐标x,y x(i,1)=0.5*(n-i+1); y(i,1) = (n-i+1)+y0; for j=2:i x(i,j)=x(i,1)+(j-1)*1; y(i,j)=y(i,1); end end mm =moviein(m); for i=1:m s=rand(1,n); %产生n个随机数 xi=x(1,1); yi=y(1,1); k=1;l=1; for j=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):n],ballnum1) axis([-2,n+2,0,y0+n+1]) %画各格子的频率 mm(i)=getframe; hold off %放弃图形保持句柄 end movie(mm,1) %播放动画一次
1701009713
1701009714
运行程序输出图形如图8-11和图8-12所示。
1701009715
1701009716
1701009717
1701009718
图8-11 运动模拟1 图8-12 运动模拟2 （2）用二项分布描述Galton钉板模型，具体实现过程如下。
1701009719
1701009720
小球自上方落下，经过n个钉子。每经过一个钉子时只有两种可能结果：向左或向右，这是一个具有两个结果（成功和失败）的随机试验E，将向右视为成功，其概率为p，向左视为失败，其概率为1-p。小球碰到一个钉子下落一格，相当于进行了一次试验E，自顶端落下，碰到n个钉子，最终落到某个格子的过程，恰好相当于将试验E重复了n次，因此一次投球过程就是一个n重贝努利试验。n重贝努利试验的成功次数X正好就是小球向右移动的次数，是一个随机变量，根据概率论的结果，它服从二项分布，即X～B（n，p）。其取值与模拟模型的对应关系如表8-4所示。
1701009721
1701009722
表8-4 格子编号与随机变量取值对应表
1701009723
1701009724
1701009725
1701009726
1701009727
利用概率论知识，二项随机变量X的分布列为：
1701009728
1701009729
1701009730
1701009731
1701009732
有了上面的理论分析之后，我们可以比较n次投球小球堆积的频率图和X～B（n，0.5）的分布图之间的差异。
1701009733
1701009734
编写MATLAB程序如下：
1701009735
1701009736
clc,clear,close all %清屏和清除变量 warning off %消除警告 %Galton钉板试验 %设置参数 m = 500; n=5; p=0.5; rand(‘seed’,3); R=binornd(n,p,1,m); %模拟服从二项分布的随机数，相当于模拟投球m次 %确定落在编号为i-1的格子的小球频率 for i=1:n+1 k=[]; k = find(R==(i-1)); h(i)=length(k)/m; end x=0:n; figure(‘color’,[1,1,1]) subplot(121) %画频率图 axis([-1,6,0,1]) bar(x,h) xlabel(‘(1)500次投球小球堆积的频率图’) f = binopdf(x,n,p); subplot(122) axis([-1,6,0,1]) bar(x,f) xlabel(‘(2)理论分布B(5,0.5)的分布图’)
1701009737
1701009738
运行程序输出图形如图8-13所示。
1701009739
1701009740
1701009741
1701009742
1701009743
图8-13 用二项分布描述Galton钉板试验
1701009744
1701009745
500次投球小球堆积频率图和理论分布X～B（n，0.5）的分布图是接近的，因此可知道，Galton试验过程中，小球运动是服从二项分布X～B（n，p）的。
1701009746
1701009747
1701009748
1701009749