概率统计的数值实验——MATLAB在概率统计教学中的应用

概率统计的数值实验——MATLAB在概率统计教学中的应用崔明涛2012年10月11日引言而MATLAB软件具有简单易学、易操作和绘图功能强等特点，利用MATLAB软件的图形可视功能将概率统计的内容用图形表示出来，通过图形让学生加深理解，以达到事半功倍的效果。概率论与数理统计知识比较抽象，逻辑性较强。因此，建议让学生结合理论和公式推导，进行数值试验和相关调查，直观地感受数学概念和理论，从而提高学生解决实际问题的信心和能力。概率论1.rand(m,n)：生成m×n的随机矩阵，每个元素都在(0,1)间，生成方式为均匀分布。2.randn(m,n)：生成m×n的随机矩阵，每个元素都在(0,1)间，生成方式为正态分布。3.randperm(m)：生成一个1~m的随机整数排列。4.perms(1:n)：生成一个1~n的全排列，共n!个。5.取整函数系列：（1）fix(x)：截尾法取整；（2）floor(x)：退一法取整（不超过x的最大整数）；（3）ceil(x)：进一法取整（=floor(x)+1）；（4）round(x)：四舍五入法取整。6.unique(a)：合并a中相同的项。7.prod(x)：向量x的所有分量元素的积。一、MATLAB常用的与随机数产生相关的函数：示例：rand(1)%生成一个(0,1)间的随机数ans=0.8147rand(2,2)%生成一个2×2阶(0,1)间的随机数矩阵ans=0.91340.09750.63240.2785randperm(5)%生成一个1~5的随机整数排列ans=41523a=[1242332];unique(a)ans=1234例1随机投掷均匀硬币，观察国徽朝上与国徽朝下的频率。解n=3000~100000000;m=0;•fori=1:n•t=randperm(2);%生成一个1~2的随机整数排列•x=t-1;%生成一个0~1的随机整数排列•y=x(1);•ify==0;•m=m+1;•end•end•p1=m/n•p2=1-p1试验次数n300050001万2万3万国徽朝上频率0.50400.50060.48790.49990.5046国徽朝下频率0.49600.49940.51210.50010.4954试验次数n5万10万100万100万1亿国徽朝上频率0.50210.49990.49990.50010.5000国徽朝下频率0.49790.50010.50010.49990.5000APAfnn可见当时，解记事件为第i个人拿到自已枪，事件为第i个人没拿到自己枪，易知：又记为没有一个人拿到自己枪的概率。有乘法公式可知：例2某班有n个人，每人各有一支枪，这些枪外形一样。某次夜间紧急集合，若每人随机地取走一支枪，问没有一个人拿到自己枪的概率是多少？iAiAnAPi1nnAPi1ni,,2,10p!1321nAAAAPnnjinnAAPAPAAPijiji111nkjinnnAAAPAPAAAPkjikkji1211于是所以特别地，当n较大时，。因此，可随机模拟出没有人拿到自己枪的频率，根据频率的稳定性，近似当做概率，然后去估计自然常数e。算法如下：!1,21\1,132131211nAAAAPnnnCAAAPnnCAAPAPnnnnkjikjinnnjijiniinkknnnniiknnnnCnnCAPp013210!1!121111110ep1、产生n个随机数的随机序列；2、检验随机列与自然列是否至少有一个配对；3、对没有一个配对的序列进行累积p；4、重复1、2、3步m次；5、估计。pme具体程序及相关结果为（注：自然常数e≈2.7183）：m=40000;n=50;p=0;forj=1:mk=0;sui=randperm(n);fori=1:nifsui(i)==ik=k+1;elsek=k;endendifk==0p=p+1;elsep=p;endende=m/pe=2.7313模拟次数m400004000040000人数n100020005000e2.71552.70822.7202模拟次数m400040000400000人数n505050e2.73792.73132.7194设针与平行线的夹角为，针的中心与最近直线的距离为。针与平行线相交的充要条件是，则所求概率为故可得的近似计算公式，其中n为随机试验次数，m为针与平行线相交的次数。例3Buffon投针实验在画有许多间距为的等距平行线的白纸上，随机投掷一根长为的均匀直针，求针与平行线相交的概率，并计算的近似值。)(dlld)0()20(dyysin2lynmdlddlp22sin20mdnl2解clear,clf•n=10000000;l=0.5;m=0;d=1;•fori=1:n•x=l/2*sin(rand(1)*pi);y=rand(1)*d/2;•ifx=y•m=m+1;•end•end•p1=m/n•pai=2*n*l/(m*d)试验次数n5千1万10万100万1000万针长l/平行间距d3/103/103/103/103/10相交频率0.18360.19710.18870.19050.1912π的近似值3.26803.04413.17983.14983.1387试验次数n5千1万10万100万1000万针长l/平行间距d2/52/52/52/52/5相交频率0.24960.25620.25490.25440.2543π的近似值3.20513.12263.13863.14513.1433试验次数n5千1万10万100万1000万针长l/平行间距d1/21/21/21/21/2相交频率0.32540.31480.31580.31780.3183π的近似值3.07313.17663.16673.14703.1417试验次数n5千1万10万100万1000万针长l/平行间距d4/54/54/54/54/5相交频率0.51420.51340.50860.50930.5093π的近似值3.11163.11653.14603.14183.1418试验次数n5千1万10万100万1000万针长l/平行间距d17/2017/2017/2017/2017/20相交频率0.54320.54520.54200.54120.5410π的近似值3.12963.11813.13663.14133.1426试验次数n5千1万10万100万1000万针长l/平行间距d9/109/109/109/109/10相交频率0.58600.57000.57560.57330.5731π的近似值3.07173.15793.12723.13953.1410例4在100个人的团体中，不考虑年龄差异，研究是否有两个以上的人生日相同。假设每人的生日在一年365天中的任意一天是等可能的，那么随机找n个人(不超过365人)。(1)求这n个人生日各不相同的概率是多少？从而求这n个人中至少有两个人生日相同这一随机事件发生的概率是多少？(2)近似计算在30名学生的一个班中至少有两个人生日相同的概率是多少？•解:(1)clear,clf•forn=1:100•p0(n)=prod(365:-1:365-n+1)/365^n;•p1(n)=1-p0(n);•end•p1=ones(1,100)-p0;•n=1:100;•plot(n,p0,n,p1,'--')•xlabel('人数'),ylabel('概率')•legend('生日各不相同的概率','至少两人生日相同的概率')•axis([0100-0.11.199]),gridon010203040506070809010000.20.40.60.81人数概率生日各不相同的概率至少两人生日相同的概率•p1(30)=0.7063,p1(60)=0.9941分析：在30名学生中至少两人生日相同的概率为70.63％。下面进行计算机仿真。随机产生30个正整数，代表一个班30名学生的生日，然后观察是否有两人以上生日相同。当30个人中有两人生日相同时，输出“1”，否则输出“0”。如此重复观察100次，计算出这一事件发生的频率。100f•(2)clear,clf•n=0;•form=1:100%做100次随机试验•y=0;•x=1+fix(365*rand(1,30));%产生30个随机数•fori=1:29%用二重循环寻找30个随机数中是否有相同数•forj=i+1:30•ifx(i)==x(j)•y=1;break;•end•end•end•n=n+y;%累计有两人生日相同的试验次数•end•f=n/m%计算频率•f=0.6900•f=0.7900•f=0.6700•f=0.7300•f=0.7500•f=0.6900•f=0.7200•f=0.6700•f=0.6800……重复观察，数据如下：例5Galton钉板模型和二项分布Galton钉板试验是由英国生物统计学家和人类学家Galton设计的。故而得名。通过模拟Calton钉板试验，观察和体会二项分布概率分布列的意义、形象地理解DeMoivre-Laplace中心极限定理。共15层小钉Ox-8-7-6-5-4-3-2-112345678小球最后落入的格数?W记小球向右落下的次数为则,X~X(15,0.5)b记小球向左落下的次数为则,Y~Y(15,0.5)b[sign()]/2WXYXY符号函数,大于0返回1,小于0返回-1,等于0返回0高尔顿(FrancisGalton,1822-1911)英国人类学家和气象学家Ox-8-7-6-5-4-3-2-112345678记则151kkX近似~)51,0(N),(2nnN共15层小钉1,1,kX(1,2,,15)k小球碰第层钉后向右落下k小球碰第层钉后向左落下k151,01522nn模拟Galton钉板试验的步骤：(1)确定钉子的位置：将钉子的横、纵坐标存储在两个矩阵X和Y中。(2)在Galton钉板试验中，小球每碰到钉子下落时都具有两种可能性，设向右的概率为p，向左的概率为q＝1-p，这里p=0.5，表示向左向右的机会是相同的。模拟过程如下：首先产生一均匀随机数u，这只需调用随机数发生器指令rand(m,n)。rand(m,n)指令：用来产生m×n个(0,1)区间中的随机数，并将这些随机数存于一个m×n矩阵中，每次调用rand(m,n)的结果都会不同。如果想保持结果一致，可与rand(‘seed’,s)配合使用，这里s是一个正整数，例如rand('seed',1),u=rand(1,6)u=0.51290.46050.35040.09500.43370.7092而且再次运行该指令时结果保持不变。除非重设种子seed的值，如rand('seed',2),u=rand(1,6)u=0.02580.92100.70080.19010.86730.4185这样结果才会产生变化。将[0,1]区间分成两段，区间[0,p)和[p,1]。如果随机数u属于[0,p)，让小球向右落下；若u属于[p,1]，让小球向左落下。将这一过程重复n次，并用直线连接小球落下时所经过的点，这样就模拟了小球从顶端随机地落人某一格子的过程。(3)模拟小球堆积的形状。输入扔球次数m(例如m＝50、100、500等等)，计算落在第i个格子的小球数在总球数m中所占的比例，这样当模拟结束时，就得到了频率用频率反映小球的堆积形状。(4)用如下动画指令制作动画：movien(n)：创建动画矩阵；制作动画矩阵数据；Getframe：拷贝动画矩阵；movie