《Matlab生成随机数》PPT课件.ppt

2019/12/5,生成随机数,教材,2019/12/5,2019/12/5,主要内容生成一元分布随机数生成多元分布随机数蒙特卡洛方法,2019/12/5,第一节生成一元分布随机数,2019/12/5,一、均匀分布随机数和标准正态分布随机数,调用格式：Y=randY=rand(n)Y=rand(m,n)Y=rand(mn)Y=rand(m,n,p,)Y=rand(mnp)Y=rand(size(A),1.rand函数,2019/12/5,在MATLAB7.7以前的版本中，rand函数还可以这样调用：rand(method,s)s=rand(method)其中method是字符串变量，它的可能取值如下表所列：,2019/12/5,调用格式：与rand函数类似,2.randn函数,2019/12/5,%设置随机数生成器的算法为MersenneTwister算法，初始种子为1rand(twister,1);%生成2行6列的随机数矩阵，其元素服从0,1上均匀分布x1=rand(2,6),【例4.1-1】设置随机数生成器的算法为MersenneTwister算法，生成均匀分布随机数矩阵,2019/12/5,二、RandStream类,MATLAB7.7及以后的版本中，依然支持rand函数的上述两种调用方式，但已经是过时的调用方式了，因为MATLAB7.7中对生成随机数作了重大调整，给出了RandStream（随机数流）类，通过调用类的构造函数并传递合适的参数可以创建类对象，然后调用类对象的rand，randn，randi，randperm方法生成随机数。,2019/12/5,RandStream类方法列表,1.创建RandStream类对象,2019/12/5,调用格式：s=RandStream(gentype)=RandStream(gentype,param1,val1,param2,val2,),2.RandStream函数的调用方法,gentype有6个可能的取值：mcg16807、mlfg6331_64、mrg32k3a、mt19937ar、shr3cong和swb2712，对应随机数生成器的6个不同算法，默认值为mt19937ar,2019/12/5,%创建一个RandStream类对象s，其随机数生成器的算法为mlfg6331_64，初始种子为10，对象s的randn方法的算法为Inversions=RandStream(mlfg6331_64,seed,10,RandnAlg,Inversion);%调用对象s的randn方法生成10行10列的随机数矩阵x，其元素服从标准正态分布x=s.randn(10)y=x(:);%将x按列拉长成一个列向量hist(y)%绘制频数直方图xlabel(标准正态分布随机数);%为X轴加标签ylabel(频数);%为Y轴加标签,【例4.1-2】调用RandStream函数创建一个指定随机数生成算法的RandStream类对象，然后利用对象的randn方法生成1010的标准正态分布随机数矩阵，并将矩阵按列拉长，画出频数直方图,2019/12/5,三、常见一元分布随机数,MATLAB统计工具箱中函数名以rnd三个字符结尾的函数用来生成常见分布的随机数。例如：betarndBeta分布exprnd指数分布gamrndGamma分布lognrnd对数正态分布normrnd正态分布poissrnd泊松分布randsample从有限总体中随机抽样random指定分布,2019/12/5,%调用normrnd函数生成1000行3列的随机数矩阵x，其元素服从均值为75，标准差为8的正态分布x=normrnd(75,8,1000,3);hist(x)%绘制矩阵x每列的频数直方图xlabel(正态分布随机数（mu=75,sigma=8）);%为X轴加标签ylabel(频数);%为Y轴加标签legend(第一列,第二列,第三列)%为图形加标注框,【例4.1-3】调用normrnd函数生成10003的正态分布随机数矩阵，其中均值为75，标准差为8，并作出各列的频数直方图,2019/12/5,%调用normrnd函数生成1000行3列的随机数矩阵x，其各列元素分别服从不同的正态分布x=normrnd(repmat(01540,1000,1),repmat(123,1000,1),1000,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，并作出各列的频数直方图,2019/12/5,%调用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轴加标签,【例4.1-5】调用random函数生成100001的二项分布随机数向量，然后作出频率直方图。其中二项分布的参数为n=10，p=0.3,2019/12/5,x=random(chi2,10,10000,1);fp,xp=ecdf(x);%计算经验累积概率分布函数值ecdfhist(fp,xp,50);%绘制频率直方图holdont=linspace(0,max(x),100);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,2019/12/5,四、任意一元分布随机数,MATLAB统计工具箱中的randsample和randsrc函数用来生成指定离散分布随机数,1.离散分布随机数,2019/12/5,xvalue=-2-1012;%定义向量xvaluexp=0.050.20.50.20.05;%定义向量xp%调用randsample函数生成100个服从指定离散分布的随机数x=randsample(xvalue,100,true,xp);reshape(x,1010)%调用randsrc函数生成10*10的服从指定离散分布的随机数矩阵y=randsrc(10,10,xvalue;xp),2019/12/5,xvalue=ABCDE;%定义向量xvaluexp=0.30.20.250.20.05;%定义向量xp%调用randsample函数生成100个服从指定离散分布的随机字符序列x=randsample(xvalue,100,true,xp);,2019/12/5,根据以上原理自编函数crnd，代码略,2.连续分布随机数,2019/12/5,pdffun=6*x*(1-x);%密度函数表达式%调用crnd函数生成1000个服从指定一元连续分布的随机数x=crnd(pdffun,01,1000,1);fp,xp=ecdf(x);%计算经验累积概率分布函数值ecdfhist(fp,xp,20);%绘制频率直方图holdonfplot(pdffun,01,r)%绘制真实密度函数曲线xlabel(x);%为X轴加标签ylabel(f(x);%为Y轴加标签legend(频率直方图,密度函数曲线)%为图形加标注框,2019/12/5,rand(twister,1);x=unifrnd(0,6,1000,1);rand(twister,1);pdffun=1/6;%密度函数表达式y=crnd(pdffun,06,1000,1);scatterhist(x,y),2019/12/5,第二节生成多元分布随机数,2019/12/5,生成多元分布随机数的MATLAB函数,2019/12/5,2019/12/5,mu=1020;%二元正态分布的均值向量sigma=13;316;%二元正态分布的协方差矩阵%调用mvnrnd函数生成10000组二元正态分布随机数xy=mvnrnd(mu,sigma,10000);hist3(xy,15,15);%绘制二元正态分布随机数的频数直方图xlabel(X)%为X轴加标签ylabel(Y)%为Y轴加标签zlabel(频数)%为Z轴加标签,2019/12/5,第三节蒙特卡洛方法,2019/12/5,蒙特卡洛(MonteCarlo)方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战期间研制原子弹的“曼哈顿计划”。该计划的主持人之一，数学家冯诺伊曼用摩纳哥的驰名世界的赌城MonteCarlo来命名这种方法，因此称之为MonteCarlo方法。,2019/12/5,一、有趣的蒙提霍尔问题,蒙提霍尔问题（MontyHallproblem），也称为三门问题，是一个源自博弈论的数学游戏问题，问题的名字来自美国的电视游戏节目：LetsMakeaDeal，该节目的主持人名叫蒙提霍尔（MontyHall）。这个游戏的玩法是：参赛者面前有三扇关闭的门，其中一扇门的后面藏有一辆汽车，而另外两扇门的后面则各藏有一只山羊。参赛者从三扇门中随机选取一扇，若选中后面有车的那扇门就可以赢得该汽车。当参赛者选定了一扇门，但尚未开启它的时候，节目主持人会从剩下两扇门中打开一扇藏有山羊的门，然后问参赛者要不要更换自己的选择，选取另一扇仍然关上的门。这个游戏涉及到的问题是：参赛者更换自己的选择是否会增加赢得汽车的概率？,2019/12/5,1.理论求解,2019/12/5,2.蒙特卡洛方法求解,functionp=SheepAndCar(n)%p=SheepAndCar(n)，利用蒙特卡洛方法求解蒙提霍尔问题，%求参赛者更换选择之后%赢得汽车的概率p.这里的n是正整数标量或向量，表示随机抽样的次数。fori=1:length(n)x=randsample(3,n(i),true);%随机抽样p(i)=sum(x=3)/n(i);%概率的模拟值end,2019/12/5,二、用蒙特卡洛方法求圆周率,1.随机投点法原理图解,2019/12/5,2.程序实现,functionpiva=PiMonteCarlo(n)%PiMonteCarlo(n)，用随机投点法模拟圆周率pi，作出模拟图.n为投点次数，可以是%非负整数标量或向量.%piva=PiMonteCarlo(n)，用随机投点法模拟圆周率pi，返回模拟值piva.若n为标量（向%量），则piva也为标量（向量）.x=0;y=0;d=0;m=length(n);%求变量n的长度pivalue=zeros(m,1);%为变量pivalue赋初值%通过循环用投点法模拟圆周率pifori=1:mx=2*rand(n(i),1)-1;y=2*rand(n(i),1)-1;d=x.2+y.2;pivalue(i)=4*sum(d1%如果n为向量，则返回圆周率的模拟值与投点个数的散点图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为标量，则返回投点法模拟圆周率的示意图figure;%新建一个图形窗口plot(x,y,k.);%绘制散点图holdon;%绘制边长为2的正方形h=rectangle(Position,-1-122,LineWidth,2);t=linspace(0,2*pi,100);%定义一个角度向量plot(cos(t),sin(t),k,linewidth,2);%绘制单位圆xlabel(X);ylabel(Y);%添加坐标轴标签title(Pi的模拟值：num2str(pivalue);%添加标题axis(-1.11.1-1.11.1);axisequal;%设置坐标轴属性endelsepiva=pivalue;%输出圆周率的模拟值end,2019/12/5,三、用蒙特卡洛方法求积分,1.定积分,2019/12/5,functionS0,Sm=quad1mont(n)%S0,Sm=quad1mont(n),求曲线y=sqrt(x)与直线y=x所围成的阴影区域的%面积的理论值S0与蒙特卡洛模拟值Sm.输入参数n是随机投点的个数，可以是正整数标%量或向量.%S0=int(sqrt(x)-x,0,1);%面积的理论值（解析解）S0=quad(x)sqrt(x)-x,0,1);%面积的理论值（数值解）%计算阴影区域的面积的蒙特卡洛模拟值fori=1:length(n)x=rand(n(i),1);%点的横坐标y=rand(n(i),1);%点的纵坐标m=sum(sqrt(x)=y%落到阴影区域内点的频率，即概率的模拟值end,2019/12/5,2.二重积分,2019/12/5,functionV0,Vm=quad2mont(n)%V0,Vm=quad2mont(n),求球面x2+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函数（matlab2009a中出现的新函数）求体积的理论值（数值解）V0=4*quad2d(x,y)sqrt(4-x.2-y.2),0,2,0,(x)sqrt(1-(1-x).2);%求体积的蒙特卡洛模拟值fori=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)%落到所求立体内的点的频率，即概率的模拟值end,2019/12/5,四、街头骗局揭秘,2019/12/5,1.理论求解,2019/12/5,2.蒙特卡洛方法求解,functionEm,E0=GameMont1(n)%摸球游戏揭秘%Em,E0=GameMont1(n)，求摸球者在一次摸球游戏中得到的奖金（罚金）的数学期%望（均值）的理论值E0和蒙特卡洛模拟值Em.输入参数n为游戏的次数，它是一个正整数.a=nchoosek(16,8);%组合数（16选8）p=0