Matlab生成随机数.ppt

《Matlab生成随机数.ppt》由会员分享，可在线阅读，更多相关《Matlab生成随机数.ppt（42页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28生成随机数生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/28教 材生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28主要内容主要内容 生成一元分布随机数生成一元分布随机数 生成多元分布随机数生成多元分布随机数 蒙特卡洛方法蒙特卡洛方法生成随机数生成随机数 谢中华谢中华,天津科技

2、大学数学系天津科技大学数学系.2023/1/282023/1/28第一节第一节 生成一元分布随机数生成一元分布随机数生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28一、均匀分布随机数和标准正态分布随机数一、均匀分布随机数和标准正态分布随机数调用格式：调用格式：Y=randY=rand(n)Y=rand(m,n)Y=rand(m n)Y=rand(m,n,p,)Y=rand(m n p)Y=rand(size(A)1.rand函数函数生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/

3、1/28在在MATLAB7.7以前的版本中，以前的版本中，rand函数还可以这样调用：函数还可以这样调用：rand(method,s)s=rand(method)其中其中method是字符串变量，它的可能取值如下表所列：是字符串变量，它的可能取值如下表所列：生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28调用格式：调用格式：与与rand函数类似函数类似2.randn函数函数生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28%设置随机数生成器的算法为设置随机数生成器的算法为M

4、ersenne Twister算法，初始种算法，初始种子为子为1 rand(twister,1);%生成生成2行行6列的随机数矩阵，其元素服从列的随机数矩阵，其元素服从0,1上均匀分布上均匀分布 x1=rand(2,6)【例例4.1-1】设置随机数生成器的算法为设置随机数生成器的算法为Mersenne Twister算法，算法，生成均匀分布随机数矩阵生成均匀分布随机数矩阵生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28二、二、RandStream类类MATLAB 7.7及以后的版本中，依然支持及以后的版本中，依然支持rand函数的函

5、数的上述两种调用方式，但已经是过时的调用方式了，因上述两种调用方式，但已经是过时的调用方式了，因为为MATLAB 7.7中对生成随机数作了重大调整，给出中对生成随机数作了重大调整，给出了了RandStream（随机数流）类，通过调用类的构造（随机数流）类，通过调用类的构造函数并传递合适的参数可以创建类对象，然后调用类函数并传递合适的参数可以创建类对象，然后调用类对象的对象的rand，randn，randi，randperm方法生成随方法生成随机数。机数。生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28RandStream类方法列表类

6、方法列表1.创建创建RandStream类对象类对象生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28调用格式：调用格式：s=RandStream(gentype)=RandStream(gentype,param1,val1,param2,val2,)2.RandStream函数的调用方法函数的调用方法gentype有有6个可能的取值：个可能的取值：mcg16807、mlfg6331_64、mrg32k3a、mt19937ar、shr3cong和和swb2712，对应随机数生，对应随机数生成器的成器的6个不同算法，默认值为个不同算法

7、，默认值为mt19937ar生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28%创建一个创建一个RandStream类对象类对象s，其随机数生成器的算法为，其随机数生成器的算法为mlfg6331_64，初始种子为，初始种子为10，对象，对象s的的randn方法的算法为方法的算法为Inversion s=RandStream(mlfg6331_64,seed,10,RandnAlg,Inversion);%调用对象调用对象s的的randn方法生成方法生成10行行10列的随机数矩阵列的随机数矩阵x，其元素服从标准正，其元素服从标准正态分布

8、态分布 x=s.randn(10)y=x(:);%将将x按列拉长成一个列向量按列拉长成一个列向量 hist(y)%绘制频数直方图绘制频数直方图 xlabel(标准正态分布随机数标准正态分布随机数);%为为X轴加标签轴加标签 ylabel(频数频数);%为为Y轴加标签轴加标签【例例4.1-2】调用调用RandStream函数创建一个指定随机数生成算法函数创建一个指定随机数生成算法的的RandStream类对象，然后利用对象的类对象，然后利用对象的randn方法生成方法生成1010的的标准正态分布随机数矩阵，并将矩阵按列拉长，画出频数直方图标准正态分布随机数矩阵，并将矩阵按列拉长，画出频数直方图生

9、成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28三、常见一元分布随机数三、常见一元分布随机数MATLAB统计工具箱中函数名以统计工具箱中函数名以rnd三个字符结尾的三个字符结尾的函数用来生成常见分布的随机数。例如：函数用来生成常见分布的随机数。例如：betarndBeta分布分布exprnd指数分布指数分布gamrndGamma分布分布lognrnd对数正态分布对数正态分布normrnd正态分布正态分布poissrnd泊松分布泊松分布randsample从有限总体中随机抽样从有限总体中随机抽样random指定分布指定分布生成随机数生成

10、随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28%调用调用normrnd函数生成函数生成1000行行3列的随机数矩阵列的随机数矩阵x，其元素服从均值为，其元素服从均值为75，标准差为，标准差为8的正态分布的正态分布 x=normrnd(75,8,1000,3);hist(x)%绘制矩阵绘制矩阵x每列的频数直方图每列的频数直方图 xlabel(正态分布随机数（正态分布随机数（mu=75,sigma=8）);%为为X轴加标签轴加标签 ylabel(频数频数);%为为Y轴加标签轴加标签 legend(第一列第一列,第二列第二列,第三列第三列)%为图形加

11、标注框为图形加标注框【例例4.1-3】调用调用normrnd函数生成函数生成10003的正态分布随机数矩阵，的正态分布随机数矩阵，其中均值其中均值 为为75，标准差为，标准差为8，并作出各列的频数直方图，并作出各列的频数直方图生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28%调用调用normrnd函数生成函数生成1000行行3列的随机数矩阵列的随机数矩阵x，其各列元素分别服从不同，其各列元素分别服从不同的正态分布的正态分布 x=normrnd(repmat(0 15 40,1000,1),repmat(1 2 3,1000,1),1

12、000,3);hist(x,50)%绘制矩阵绘制矩阵x每列的频数直方图每列的频数直方图 xlabel(正态分布随机数正态分布随机数);%为为X轴加标签轴加标签 ylabel(频数频数);%为为Y轴加标签轴加标签%为图形加标注框为图形加标注框 legend(mu=0,sigma=1,mu=15,sigma=2,mu=40,sigma=3)【例例4.1-4】调用调用normrnd函数生成函数生成10003的正态分布随机数矩阵，的正态分布随机数矩阵，其中第各列均值分别为其中第各列均值分别为0，15，40，标准差分别为，标准差分别为1，2，3，并作，并作出各列的频数直方图出各列的频数直方图生成随机数生

13、成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28%调用调用random函数生成函数生成10000行行1列的随机数向量列的随机数向量x，其元素服从二项分布，其元素服从二项分布B(10,0.3)x=random(bino,10,0.3,10000,1);fp,xp=ecdf(x);%计算经验累积概率分布函数值计算经验累积概率分布函数值 ecdfhist(fp,xp,50);%绘制频率直方图绘制频率直方图 xlabel(二项分布（二项分布（n=10,p=0.3）随机数）随机数);%为为X轴加标签轴加标签 ylabel(f(x);%为为Y轴加标签轴加标

14、签【例例4.1-5】调用调用random函数生成函数生成100001的二项分布随机数向量，的二项分布随机数向量，然后作出频率直方图。其中二项分布的参数为然后作出频率直方图。其中二项分布的参数为 n=10，p=0.3 生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28 x=random(chi2,10,10000,1);fp,xp=ecdf(x);%计算经验累积概率分布函数值计算经验累积概率分布函数值 ecdfhist(fp,xp,50);%绘制频率直方图绘制频率直方图 hold on t=linspace(0,max(x),100);

15、y=chi2pdf(t,10);plot(t,y,r,linewidth,3)xlabel(x (chi2(10);%为为X轴加标签轴加标签 ylabel(f(x);%为为Y轴加标签轴加标签 legend(频率直方图频率直方图,密度函数曲线密度函数曲线)%为图形加标注框为图形加标注框【例例4.1-6】调用调用random函数生成函数生成100001的卡方分布随机数向量，的卡方分布随机数向量，然后作出频率直方图，并与自由度为然后作出频率直方图，并与自由度为10的卡方分布的密度函数曲的卡方分布的密度函数曲线作比较。其中卡方分布的参数（自由度）为线作比较。其中卡方分布的参数（自由度）为 10生成随机

16、数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28四、任意一元分布随机数四、任意一元分布随机数MATLAB统计工具箱中的统计工具箱中的randsample和和randsrc函数函数用来生成指定离散分布随机数用来生成指定离散分布随机数1.离散分布随机数离散分布随机数生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28 xvalue=-2-1 0 1 2;%定义向量定义向量xvalue xp=0.05 0.2 0.5 0.2 0.05;%定义向量定义向量xp%调用调用randsample函

17、数生成函数生成100个服从指定离散分布的随机数个服从指定离散分布的随机数 x=randsample(xvalue,100,true,xp);reshape(x,10 10)%调用调用randsrc函数生成函数生成10*10的服从指定离散分布的随机数矩阵的服从指定离散分布的随机数矩阵 y=randsrc(10,10,xvalue;xp)生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28 xvalue=ABCDE;%定义向量定义向量xvalue xp=0.3 0.2 0.25 0.2 0.05;%定义向量定义向量xp%调用调用randsa

18、mple函数生成函数生成100个服从指定离散分布的随机字符序列个服从指定离散分布的随机字符序列 x=randsample(xvalue,100,true,xp);生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28根据以上原理自编函数根据以上原理自编函数crnd，代码略，代码略2.连续分布随机数连续分布随机数生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28 pdffun=6*x*(1-x);%密度函数表达式密度函数表达式%调用调用crnd函数生成函数生成1000个服从指定一元

19、连续分布的随机数个服从指定一元连续分布的随机数 x=crnd(pdffun,0 1,1000,1);fp,xp=ecdf(x);%计算经验累积概率分布函数值计算经验累积概率分布函数值 ecdfhist(fp,xp,20);%绘制频率直方图绘制频率直方图 hold on fplot(pdffun,0 1,r)%绘制真实密度函数曲线绘制真实密度函数曲线 xlabel(x);%为为X轴加标签轴加标签 ylabel(f(x);%为为Y轴加标签轴加标签 legend(频率直方图频率直方图,密度函数曲线密度函数曲线)%为图形加标注框为图形加标注框生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科

20、技大学数学系.2023/1/282023/1/28 rand(twister,1);x=unifrnd(0,6,1000,1);rand(twister,1);pdffun=1/6;%密度函数表达式密度函数表达式 y=crnd(pdffun,0 6,1000,1);scatterhist(x,y)生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28第二节第二节 生成多元分布随机数生成多元分布随机数生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28生成多元分布随机数的生成多元分布

21、随机数的MATLAB函数函数生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28 mu=10 20;%二元正态分布的均值向量二元正态分布的均值向量 sigma=1 3;3 16;%二元正态分布的协方差矩阵二元正态分布的协方差矩阵%调用调用mvnrnd函数生成函数生成10000组二元正态分布随机数组二元正态分布随机数 xy=mvnrnd(mu,sigma,10000);hist3(xy,15,15);%绘制二元正态分布随机数的频数

22、直方图绘制二元正态分布随机数的频数直方图 xlabel(X)%为为X轴加标签轴加标签 ylabel(Y)%为为Y轴加标签轴加标签 zlabel(频数频数)%为为Z轴加标签轴加标签生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28第三节第三节 蒙特卡洛方法蒙特卡洛方法生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28蒙特卡洛蒙特卡洛(Monte Carlo)方法，或称计方法，或称计算机随机模拟方法，是一种基于算机随机模拟方法，是一种基于“随随机数机数”的计算方法。这一方法源于美

23、的计算方法。这一方法源于美国在第一次世界大战期间研制原子弹国在第一次世界大战期间研制原子弹的的“曼哈顿计划曼哈顿计划”。该计划的主持人。该计划的主持人之一，数学家冯之一，数学家冯诺伊曼用摩纳哥的驰诺伊曼用摩纳哥的驰名世界的赌城名世界的赌城Monte Carlo来命名这种来命名这种方法，因此称之为方法，因此称之为Monte Carlo方法。方法。生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28一、有趣的蒙提霍尔问题一、有趣的蒙提霍尔问题蒙提霍尔问题蒙提霍尔问题（Monty Hall problem），也称为三门问题，是），也称为三门问

24、题，是一个源自博弈论的数学游戏问题，问题的名字来自美国的电视一个源自博弈论的数学游戏问题，问题的名字来自美国的电视游戏节目：游戏节目：Lets Make a Deal，该节目的主持人名叫蒙提，该节目的主持人名叫蒙提霍尔霍尔（Monty Hall）。）。这个这个游戏的玩法游戏的玩法是：参赛者面前有三扇关闭的门，其中一扇门是：参赛者面前有三扇关闭的门，其中一扇门的后面藏有一辆汽车，而另外两扇门的后面则各藏有一只山羊。的后面藏有一辆汽车，而另外两扇门的后面则各藏有一只山羊。参赛者从三扇门中随机选取一扇，若选中后面有车的那扇门就参赛者从三扇门中随机选取一扇，若选中后面有车的那扇门就可以赢得该汽车。当参

25、赛者选定了一扇门，但尚未开启它的时可以赢得该汽车。当参赛者选定了一扇门，但尚未开启它的时候，节目主持人会从剩下两扇门中打开一扇藏有山羊的门，然候，节目主持人会从剩下两扇门中打开一扇藏有山羊的门，然后问参赛者要不要更换自己的选择，选取另一扇仍然关上的门。后问参赛者要不要更换自己的选择，选取另一扇仍然关上的门。这个游戏涉及到的问题是：参赛者更换自己的选择是否会增加这个游戏涉及到的问题是：参赛者更换自己的选择是否会增加赢得汽车的概率？赢得汽车的概率？生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/281.理论求解理论求解生成随机数生成随机数

26、谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/282.蒙特卡洛方法求解蒙特卡洛方法求解function p=SheepAndCar(n)%p=SheepAndCar(n)，利用蒙特卡洛方法求解蒙提霍尔问题，利用蒙特卡洛方法求解蒙提霍尔问题，%求参赛者更换选择之后求参赛者更换选择之后%赢得汽车的概率赢得汽车的概率p.这里的这里的n是正整数标量或向量，表示随机抽样的次数。是正整数标量或向量，表示随机抽样的次数。for i=1:length(n)x=randsample(3,n(i),true);%随机抽样随机抽样 p(i)=sum(x=3)/n(i);%概率的

27、模拟值概率的模拟值end生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28二、用蒙特卡洛方法求圆周率二、用蒙特卡洛方法求圆周率1.随机投点法原理图解随机投点法原理图解生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/282.程序实现程序实现function piva=PiMonteCarlo(n)%PiMonteCarlo(n)，用随机投点法模拟圆周率，用随机投点法模拟圆周率pi，作出模拟图，作出模拟图.n为投点次数，可以是为投点次数，可以是%非负整数标量或向量非负整数标量或向量

28、.%piva=PiMonteCarlo(n)，用随机投点法模拟圆周率，用随机投点法模拟圆周率pi，返回模拟值，返回模拟值piva.若若n为标量（向为标量（向%量），则量），则piva也为标量（向量）也为标量（向量）.x=0;y=0;d=0;m=length(n);%求变量求变量n的长度的长度pivalue=zeros(m,1);%为变量为变量pivalue赋初值赋初值%通过循环用投点法模拟圆周率通过循环用投点法模拟圆周率pifor i=1:m x=2*rand(n(i),1)-1;y=2*rand(n(i),1)-1;d=x.2+y.2;pivalue(i)=4*sum(d 1%如果如果n为向

29、量，则返回圆周率的模拟值与投点个数的散点图为向量，则返回圆周率的模拟值与投点个数的散点图 figure;%新建一个图形窗口新建一个图形窗口 plot(n,pivalue,k.);%绘制散点图绘制散点图 h=refline(0,pi);%添加参考线添加参考线 set(h,linewidth,2,color,k);%设置参考线属性设置参考线属性 text(1.05*n(end),pi,pi,fontsize,15);%添加文本信息添加文本信息 xlabel(投点个数投点个数);ylabel(pi的模拟值的模拟值);%添加坐标轴标签添加坐标轴标签 else%如果如果n为标量，则返回投点法模拟圆周率的

30、示意图为标量，则返回投点法模拟圆周率的示意图 figure;%新建一个图形窗口新建一个图形窗口 plot(x,y,k.);%绘制散点图绘制散点图 hold on;%绘制边长为绘制边长为2的正方形的正方形 h=rectangle(Position,-1-1 2 2,LineWidth,2);t=linspace(0,2*pi,100);%定义一个角度向量定义一个角度向量 plot(cos(t),sin(t),k,linewidth,2);%绘制单位圆绘制单位圆 xlabel(X);ylabel(Y);%添加坐标轴标签添加坐标轴标签 title(Pi的模拟值：的模拟值：num2str(pivalu

31、e);%添加标题添加标题 axis(-1.1 1.1-1.1 1.1);axis equal;%设置坐标轴属性设置坐标轴属性 endelse piva=pivalue;%输出圆周率的模拟值输出圆周率的模拟值end生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28三、用蒙特卡洛方法求积分三、用蒙特卡洛方法求积分1.定积分定积分生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28function S0,Sm=quad1mont(n)%S0,Sm=quad1mont(n),求曲线求曲

32、线 y=sqrt(x)与直线与直线 y=x 所围成的阴影区域的所围成的阴影区域的%面积的理论值面积的理论值S0与蒙特卡洛模拟值与蒙特卡洛模拟值Sm.输入参数输入参数n是随机投点的个数，可以是正整数标是随机投点的个数，可以是正整数标%量或向量量或向量.%S0=int(sqrt(x)-x,0,1);%面积的理论值（解析解）面积的理论值（解析解）S0=quad(x)sqrt(x)-x,0,1);%面积的理论值（数值解）面积的理论值（数值解）%计算阴影区域的面积的蒙特卡洛模拟值计算阴影区域的面积的蒙特卡洛模拟值for i=1:length(n)x=rand(n(i),1);%点的横坐标点的横坐标 y=

33、rand(n(i),1);%点的纵坐标点的纵坐标 m=sum(sqrt(x)=y&y=x);%落到阴影区域内点的频数落到阴影区域内点的频数 Sm(i)=m/n(i);%落到阴影区域内点的频率，即概率的模拟值落到阴影区域内点的频率，即概率的模拟值end生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/282.二重积分二重积分生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28function V0,Vm=quad2mont(n)%V0,Vm=quad2mont(n),求球面求球面 x

34、2+y2+z2=4 被圆柱面被圆柱面x2+y2=2*x所截得所截得%的（含在圆柱面内的部分）立体的体积的理论值的（含在圆柱面内的部分）立体的体积的理论值V0与蒙特卡洛模拟值与蒙特卡洛模拟值Vm.输入参数输入参数n是是%随机投点的个数，可以是正整数标量或向量随机投点的个数，可以是正整数标量或向量.%V0=32*(pi/2-2/3)/3;%体积的理论值体积的理论值%V0=4*quadl(x)arrayfun(xx)quadl(y)sqrt(4-xx.2-y.2),.%0,sqrt(1-(1-xx).2),x),0,2);%体积的理论值（数值解）体积的理论值（数值解）%调用调用quad2d函数（函数

35、（matlab2009a中出现的新函数）求体积的理论值（数值解）中出现的新函数）求体积的理论值（数值解）V0=4*quad2d(x,y)sqrt(4-x.2-y.2),0,2,0,(x)sqrt(1-(1-x).2);%求体积的蒙特卡洛模拟值求体积的蒙特卡洛模拟值for i=1:length(n)x=2*rand(n(i),1);%点的点的x坐标坐标 y=rand(n(i),1);%点的点的y坐标坐标 z=2*rand(n(i),1);%点的点的z坐标坐标%落到区域落到区域 T 内的点的频数内的点的频数 m=sum(x.2+y.2+z.2=4)&(x-1).2+y.2=1);Vm(i)=16*

36、m/n(i);%落到所求立体内的点的频率，即概率的模拟值落到所求立体内的点的频率，即概率的模拟值end生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/28四、街头骗局揭秘四、街头骗局揭秘生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/281.理论求解理论求解生成随机数生成随机数 谢中华谢中华,天津科技大学数学系天津科技大学数学系.2023/1/282023/1/282.蒙特卡洛方法求解蒙特卡洛方法求解function Em,E0=GameMont1(n)%摸球游戏揭秘摸球游戏揭秘

37、%Em,E0=GameMont1(n)，求摸球者在一次摸球游戏中得到的奖金（罚金）的数学期，求摸球者在一次摸球游戏中得到的奖金（罚金）的数学期%望（均值）的理论值望（均值）的理论值E0和蒙特卡洛模拟值和蒙特卡洛模拟值Em.输入参数输入参数n为游戏的次数，它是一个正整数为游戏的次数，它是一个正整数.a=nchoosek(16,8);%组合数（组合数（16选选8）p=0;%通过循环计算每种可能结果的理论概率通过循环计算每种可能结果的理论概率for i=4:8 p(i-3)=2(i=4)*nchoosek(8,i)*nchoosek(8,8-i)/a;endE0=p*-3,0.2,0.5,1,10;

38、%数学期望的理论值数学期望的理论值Freq0=zeros(1,5);%做做n次摸球游戏，计算各种可能的结果出现的频数次摸球游戏，计算各种可能的结果出现的频数for i=1:n x=randsample(16,8,false);%不放回随机抽样不放回随机抽样 x(x8)=2;%将将x中取值为中取值为9至至16的元素改为的元素改为2,用来标记白球用来标记白球%统计统计x中中1和和2出现的次数，整理成出现的次数，整理成4 4,3 5,2 6,1 7,0 8的形式的形式 x=sort(x);%将将x从小到大排序从小到大排序 x1=diff(x);%求差分求差分 x1(end+1)=1;%在在x1的最后

39、补上一个的最后补上一个1 x1=find(x1);%查找查找x1中非零元素的下标中非零元素的下标 x1=0;x1;%在在x1的前面补上一个的前面补上一个0 Freq=sort(diff(x1);%对对x1求差分，然后排序求差分，然后排序%统计统计4 4,3 5,2 6,1 7,0 8各自出现的频数各自出现的频数 if Freq=8 Freq=0;8;end if isequal(Freq,4;4)Freq0=Freq0+1 0 0 0 0;%计算计算 4:4 出现的频数出现的频数 elseif isequal(Freq,3;5)Freq0=Freq0+0 1 0 0 0;%计算计算 3:5 出现的频数出现的频数 elseif isequal(Freq,2;6)Freq0=Freq0+0 0 1 0 0;%计算计算 2:6 出现的频数出现的频数 elseif isequal(Freq,1;7)Freq0=Freq0+0 0 0 1 0;%计算计算 1:7 出现的频数出现的频数 else Freq0=Freq0+0 0 0 0 1;%计算计算 0:8 出现的频数出现的频数 endendEm=Freq0*-3,0.2,0.5,1,10/n;%计算数学期望的模拟值计算数学期望的模拟值