Matlab高尔顿板仿真模拟实验

1.前言

学习了概率论，决定动手编程实现高尔顿板仿真模拟实验,其中本文代码部分参考了这篇文章.

2.实验目的

通过模拟实验实现高尔顿板模型，并对结果与二项分布进行比较.

3.实验原理

在高尔顿钉板中倒入小球，下落过程中小球掉在左边和掉在右边的概率相同.放入若干小球，当小球数量足够多时，落下的小球所呈现的曲线符合二项分布概率.

4.实验过程

本实验采用动画的形式直观的表现实验过程，并通过柱状图表现每一次实验结果，方便与二项分布曲线进行比较，实验设置高尔顿板为7层，函数将通过生成随机数来模拟每一次小球下落的结果，并通过更改小球个数$n$来得到不同的结果.

5.实验结果

$n=100$时:

$n=1000$时:

$n=10000$时:

6.实验结论

对结果进行比较可知，n取100时初步体现了二项分布的特征，n取1000时结果更加接近二项分布曲线，当n取到10000时实验结果几乎与二项分布曲线重合.对此，我们可以推测，当试验次数足够多时，即n足够大时，高尔顿版实验结果满足二项分布规律。

7.Matlab代码

K=10000;  %控制实验的次数，可更改
n=zeros(1,7);  % 矩阵n为1*5的矩阵初值为0
delay=0.0005;  %设置delay变量等于0.005
clf;  %清除当前窗口的图案
axis([-8 8 -15 0]);  %设置x范围-4到4 y范围-7到0
aa=sqrt(3)/2;  %aa=（根号3）/2
bb=1/2;  %bb= 1/2
dd=3/2;  %d= 3/2
a=0.9*aa;  %a= 0.9 * aa
b=0.9*bb;  %b= 0.9 * bb
xxx=[-3*sqrt(3),-2*sqrt(3),-1*sqrt(3),0,1*sqrt(3),2*sqrt(3),3*sqrt(3)]; % 设定柱形图横坐标
xx=[0,-aa,aa,-2*aa,0,2*aa,-3*aa,-aa,aa,3*aa,-4*aa,-2*aa,0,2*aa,4*aa,-5*aa,-3*aa,-aa,aa,3*aa,5*aa];  % 每个六边形顶点的坐标
yy=[0,-dd,-dd,-2*dd,-2*dd,-2*dd,-3*dd,-3*dd,-3*dd,-3*dd,-4*dd,-4*dd,-4*dd,-4*dd,-4*dd,-5*dd,-5*dd,-5*dd,-5*dd,-5*dd,-5*dd];  % 每个六边形顶点的坐标
for k=1:K  %循环K次
   figure(1);
   clf  %清除当前窗口的图案
   r=round(rand(1,6));  %round函数四舍五入取整 rand函数生成一个1行4列的随机数
   x0=0;y0=0;    %
   for i=1:21  % i 循环 从1 到 10,10个六边形
      x=xx(i);  % x等于每一个六边形的顶端
      y=yy(i)-0.1;  %y等于每一个六边形的顶端
      X=[x,x-a,x-a,x,x+a,x+a];  % 标出六边形的六个点的坐标
      Y=[y,y-b,y-b-0.9,y-1.8,y-b-0.9,y-b];  % 标出六边形的六个点的坐标
      patch(X,Y,'g');hold on;  %用绿色把XY中的部分填充
   end  
   n(sum(r)+1)=n(sum(r)+1)+1;  %
   plot(x0,y0,'mo');  %画（0，0）的红圈圈
   for j=1:6  % j 循环 从1到4
      if r(j)==0  %如果r矩阵中的第j个元素为0
         x0=x0-aa;y0=y0-bb;  % 球往左移动
         plot(x0,y0,'mo');hold on;  %画红圈圈，红圈圈暂时保留，暂停0.0005s
         y0=y0-1;  %y0往下一单位距离，即往下走的圈圈
         plot(x0,y0,'mo');hold on;  %画红圈圈，红圈圈暂时保留，暂停0.0005s
      else  %如果r矩阵中的第j个元素为1
         x0=x0+aa;y0=y0-bb;  % 球往右移动
         plot(x0,y0,'mo');hold on;  %画红圈圈，红圈圈暂时保留，暂停0.0005s
         y0=y0-1;  %y0往下一单位距离，即往下走的圈圈
         plot(x0,y0,'mo');hold on;  %画红圈圈，红圈圈暂时保留，暂停0.0005s
      end  
   end  
   for m=1:7  % m 循环 从1到5
      text((m-4)*sqrt(3),-9.5,num2str(n(m)));  % 最底端的位置显示字符串（掉落球的实时结果）
   end
   figure(2);
   bar(xxx,n); % 柱形图
   figure(1);  
   pause(1*delay);  %暂停 初始时间20*0.0005=0.001s
end