matlab第m列第n个值,第12章-MATLAB-数值模拟实例解析.ppt

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

【例12-3】设袋中有10个球，其中3个是红球，7个是白球，从中取5个球，设X表示取到红球的个数，求X的概率分布。 解：易知X服从参数为10,3,5的超几何分布，故编写如下简单程序： for k=1:4 p(k)=hygepdf(k-1,10,3,5); % 生成服从参数为10，3，5的超几何分布随机数 end p 下面利用蒙特卡罗方法求解该例，编写如下语句： rand('state',2009) % 设置随机数状态 s1=0;s2=0;s3=0;s4=0; % 设置计数器 N=50000; % 模拟次数 for i=1:N x=randperm(10); % 产生1：10的一个排列 if sum(x(1:5)<=3)==0 % 取到0个红球的情形 s1=s1+1; elseif sum(x(1:5)<=3)==1 % 取到1个红球的情形 s2=s2+1; elseif sum(x(1:5)<=3)==2 % 取到两个红球的情形 s3=s3+1; else % 取到3个红球的情形 s4=s4+1; end end p=[s1 s2 s3 s4]/N 【例12-2】 准则演示。 解：以为例绘制图形进行说明，具体程序代码如下： mu=3;sigma=0.5; % 正态分布参数设定 x=mu+sigma*[-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') % 绘制图形 hold on plot([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)],[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))) 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-3\sigma','\mu-2\sigma','\mu-\sigma','\mu',... '\mu+\sigma','\mu+2\sigma','\mu+3\sigma'}) title('3\sigma准则图形表示') 【例12-5】概率分布交互界面演示。 解：在MATLAB命令窗口中输入disttool。 【例12-8】用蒙特卡罗法计算定积分和重积分。 解：①编写如下语句： s=[];a=-1;b=1;N=100000; for k=1:6 x=rand(N,2); s=[s sum(exp(-((-1+2*x(:,1))).^2)-x(:,2)>=0)/N]; end s=s*(b-a)*1 ② N=1000000; for k=1:6 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);