第12章 MATLAB 数值模拟实例解析课件.ppt

《第12章 MATLAB 数值模拟实例解析课件.ppt》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、实例解析【例12-3】设袋中有10个球，其中3个是红球，7个是白球，从中取5个球，设X表示取到红球的个数，求X的概率分布。解：易知X服从参数为10,3,5的超几何分布，故编写如下简单程序：fork=1:4p(k)=hygepdf(k-1,10,3,5);%生成服从参数为10，3，5的超几何分布随机数endp运行结果：p=0.08330.41670.41670.0833下面利用蒙特卡罗方法求解该例，编写如下语句：rand('state',2009)%设置随机数状态s1=0;s2=0;s3=0;s4=0;%设置计数器N=50000;%模拟次数fori=1:Nx=

2、randperm(10);%产生1：10的一个排列ifsum(x(1:5)<=3)==0%取到0个红球的情形s1=s1+1;elseifsum(x(1:5)<=3)==1%取到1个红球的情形s2=s2+1;elseifsum(x(1:5)<=3)==2%取到两个红球的情形s3=s3+1;else%取到3个红球的情形s4=s4+1;endendp=[s1s2s3s4]/N运行结果：p=0.08230.41810.41730.0824【例12-2】准则演示。解：以为例绘制图形进行说明，具体程序代码如下：mu=3;sigma=0.5;%正态分布参数设定x=mu+s

3、igma*[-3:-1,1:3];yf=normcdf(x,mu,sigma);%计算累积分布函数值yh=normpdf(x,mu,sigma);%计算概率密度值P=[yf(4)-yf(3),yf(5)-yf(2),yf(6)-yf(1)];xd=1:0.1:5;yd=normpdf(xd,mu,sigma);%计算概率密度值plot(xd,yd,'b')%绘制图形holdonplot([x(1),x(1)],[0,yh(1)],[x(6),x(6)],[0,yh(6)],...[x(1),2.7],[yh(1),yh(1)],'k:',[3.15,x(6)

4、],[yh(1),yh(1)],'k:')plot([x(2),x(2)],[0,yh(2)],[x(5),x(5)],[0,yh(5)],...[x(2),2.7],[yh(2),yh(2)],'k:',[3.15,x(5)],[yh(2),yh(2)],'k:')plot([x(3),x(3)],[0,yh(3)],[x(4),x(4)],[0,yh(4)],...[x(3),2.7],[yh(3),yh(3)],'k:',[3.15,x(4)],[yh(3),yh(3)],'k:')text(mu-0.5*sigma,yh(3),num2str(P(1

5、)))text(mu-0.5*sigma,yh(5),num2str(P(2)))text(mu-0.5*sigma,yh(1),num2str(P(3)))set(gca,'xticklabel',[],'yticklabel',[])x=[x(1)-0.1,x(2)-0.1,x(3)-0.1,mu,x(4)-0.1,x(5)-0.1,x(6)-0.1];y=-0.02*ones(1,7);text(x,y,{'mu-3sigma','mu-2sigma','mu-sigma','mu',...'mu+sigma','mu+2sigm

6、a','mu+3sigma'})title('3sigma准则图形表示')【例12-5】概率分布交互界面演示。解：在MATLAB命令窗口中输入disttool。【例12-8】用蒙特卡罗法计算定积分和重积分。解：①编写如下语句：s=[];a=-1;b=1;N=100000;fork=1:6x=rand(N,2);s=[ssum(exp(-((-1+2*x(:,1))).^2)-x(:,2)>=0)/N];ends=s*(b-a)*1运行结果：s=1.49311.49151.49211.49521.49351.4947②N=1000000;fork=1:6

7、x=unifrnd(-1,1,N,1);%产生区间[-1,1]上的N个均匀随机数y=unifrnd(-1,1,N,1);%产生区间[-1,1]上的N个均匀随机数z=unifrnd(0,2,N,1);%产生区间[0,2]上的N个均匀随机数Z=1+sqrt(1-x.^2-y.^2)-sqrt(x.^2+y.^2);p(k)=8*sum(z<=Z)/N;endp运行结果：p=3.14323.13993.14273.14703.14093.1480③利用有利随机数法求解该题。首先估算它的近似解：因此取有利随机数密度函数根据前面的介绍编写如下语句：rand('stat

8、e',20)%设置随机数状态fork=1:6j=1;