概率论问题MATLAB仿真求解程序

x=unifrnd(0,d/2); %随机产生数的长度，即投针之后针中点与平行线的距离 p=unifrnd(0,pi); if x<=L*sin(p)/2 m=m+1; else end end p=m/n; P(1,i)=p; i=i+1; end P=P; Pi_m=Pi_m; P_mean=mean(P) %n 次中与平行线相交的次数的频率比，即相交的概率 %记录第 i 次循环之后的相交概率值 %进入下次循环迭代 %无“;”则显示每次的相交概率值 %无“;”则显示每次的圆周率 pi 值 %显示 N 次迭代之后的相交概率均值 pi_m=(2*L*n)/(m*d); %利用投针频率估计圆周率 pi Pi_m(1,i)=(pi_m); %记录第 i 次循环之后的圆周率 pi 值
频数直方图模拟源程序
clc; close all; clear; rand('state',0) %固定随机数产生状态 R=unidrnd(2,8,10000)-1;%产生10000个1~7的随机数 test=sum(R); %对随机数求和 h=hist(test,9);%统计各数出现的频数 bpdf=binopdf(0:8,8,0.5);%计算二项分布出现各取值（即第1至第9槽内）的理论值 bpdf=bpdf/max(bpdf)*max(h);%使模拟数据与理论数据保持量级一致 b1=bar(1:9,[h;bpdf]');%绘制直方图 set(b1(2),'facecolor','r');%设置右侧的直线条填充颜色为红色 set(b1(1),'facecolor',[0.4,0.4,0.4]);%把左侧直第填充颜色设置为灰色
%n 次中与平行线相交的次数的频率比，即相交的概率,vpa() 以任意精度（4 位小数点，默认值为 32 位）显示出来 pi_m=vpa((2*L*n)/(m*d),15) %利用投针频率估计圆周率 pi,vpa()以任意精度（15 位小数点，默认值为 32 位）显示出来
运行结果
p =.3207 pi_m =3.11817898347365
k
二项分布的概率密度，可用下列程序获得：bpdf=binopdf(0:8,8,0.5) 下：
运行如
源程序 clear all; bpdf=binopdf(0:8,8,0.5) 运行结果
bpdf =0.0039 0.0312 0.1094 0.2187 0.2734 0.2187 0.1094 0.0312 0.0039 注：第1到第9个槽内的概率值，即落入第9个槽内的概率为0.0039.
（B）图源程序
%Galton板实验 clear; clc; close all; figure; xlim([-1,9]); ylim([0,9]); axis equal; hold on; L=0.8; No=8;% level number for N=1:No; mN=No/2-[1+N]/2; for k=1:N; plot([mN+k]*[1,1],[0,0.7]+No-N,'k','linewidth',2); text(-0.9,No-N+0.3,num2str(N),'fontsize',14,... 'fontname','times new roman'); end end arg=linspace(0,pi*2,200); col=[0.4,0.4,0.4]; fill(No/2+1.5+0.2*cos(arg),No/2-0.5+0.2*sin(arg),col,... 'Edgecolor',col); plot([No/2-1,No/2-1,-0.5,-0.5],... [No-1+0.7,No-1,0.7,0],'k'); plot([No/2+1,No/2+1,No+0.5,No+0.5],... [No-1+0.7,No-1,0.7,0],'k'); plot([-1,No+1],[0,0],'k','linewidth',2); for k=1:No+1; text(k-1.1,0.4,num2str(k),'fontsize',14,... 'fontname','times new roman'); plot(k-1+0.24*cos(arg),0.4+0.24*sin(arg),'k'); end set(gcf,'color','w'); axis off;
rand('state',m); %设置随机数状态 for k=1:N; at=a;%初始化甲的赌本 bt=b;%初始化乙的赌本 while at>0.5&bt>0.5;%模拟整个赌博过程 r=[(rand<p)-0.5]*2; % 算输赢 at=at+r;%交换赌本 bt=bt-r;%交换赌本 end S=S+(at<0.5); %如甲输，累加甲输的次娄 end P=S/N %计算甲输的概率值 g=p/[1-p]; Po=[1-g^b]/[1-g^(a+b)] %返回甲输光的概率理论值
MATLAB实现Buffon问题仿真求解程序
程序1
clear all; L=1; d=2; m=0; n=10000; for k=1:n x=unifrnd(0,d/2); p=unifrnd(0,pi); if x<=L*sin(p)/2 m=m+1; else end end p=vpa(m/n,4) %针的长度； %平行线间的距离（d>L）； %统计满足针与线相交条件的次数并赋初值； %投针试验次数 %迭代次数 %随机产生数的长度，即投针之后针中点与平行线的距离 %随机产生的针与线相交的角度 %针与线相交的条件 %针与线相交则记数
P =0.0676 Po =0.0656
Binomial（二项分布）的使用
Galton板实验 一个8级Galton板实验系统如下图(A)所示 (源程序见后)和(B)所示(源程序见后)
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 9
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 9
通过理论分析（见胡尧编写讲义教材）可知甲输光的概率是：wk.baidu.com
⎧ b ⎪a + b ⎪ P=⎨ 1 − ( p q )b ⎪ ⎪1 − ( p q ) a +b ⎩ 1 p=q= 2 1 p, q ≠ 2
p ,说时赌徒甲获胜,相 应甲得到一元钱,而乙付出一元钱；反之甲拿出一元钱给乙.这里对甲的赌本 a 、乙的赌 本 b 、甲赢的概率 p 取不同的数值进行 10000 次赌博过程模拟，相应程序如下：
（A）
（B）
在图（A）中，当小球从顶部向下降落时，遇到第一层竖隔板，此时小球分别向左 右下落的概率各占一半（0.5）；当小球继续下落遇到第二层竖隔板时，小球仍以 左右相同的概率往下落，以后每层均如此（如图（B）所示）。最后到了第 8 层底 部，小球将落入底部 9 个槽中的一个。但是小球落入每个槽内的概率是不一样的。 如查将这个 9 个槽编号为：1、2、3、45、6、7、8、9，计算小球落入第 9 个槽中 的概率。 理论分析：这是一个经典的二项分布概率模型。考虑在第 K 层小球运动方向有两种 可能，用 X k 表示，则 X k 服从两点分布。这里用 X k = 1 表示向右侧竖隔板方向运 动，用 X k = 0 表示向左侧竖隔板方向运动，它们发生的概率均为 0.5。最终位置 X 由 X = ∑ X k 决定，即二项分布决定，上述第 8 层即有 X ~ Binomial (8, 0.5) ，上述该
程序2
clear all; N=10; P=zeros(1,N); %循环迭代次数 %赋初值，每次循环迭之后的针与线相交的概率 p 的记录值
Pi_m=zeros(1,N); %赋初值，每次循环迭之后的圆周率 pi_m 的记录值 for i=1:N L=1; d=2; m=0; n=10000; for k=1:n %针的长度 %平行线间的距离（d>L） %统计满足针与线相交条件的次数并赋初值 %投针试验次数 %迭代次数 %随机产生的针与线相交的角度 %针与线相交的充要条件
下面用Matlab模拟小球下落到底部9个槽内的概率分布。在各层 中用0和1分别表示向左和向右运动。在小球下落到底层的槽内 时，一个8位的二进制数就完全确定了。而所落入槽的编号应该 是这个8位二进制数各位数按十进制数相加的结果。重复模拟小 球下落过程10000次，可以统计出小球落入各槽内的次数（即频 数）。画出9个频数数据的直方图。如下图所示（源程序见后）
%针与线相交则记数
运行结果
Pi_m_mean=mean(Pi_m)%显示 N 次迭代之后的圆周率 pi 均值
P_mean =0.318250000000000 Pi_m_mean =3.142648986529731
赌徒输光问题
两个赌徒甲、乙两人将进行一系列赌博。在每一局中甲获胜的概率为 p ， 而乙获胜的概率为 q ( p + q = 1 )。在每一局后，失败者都要支付一元线给 胜利者。在开始时甲拥有赌本 a 元，而乙拥有赌本 b 元，两个赌徒直到甲 输光或乙输光为止。求甲输光的概率。
3000
2500
2000
1500
1000
500
0
1
2
3
4
5
6
7
8
9
注：灰色线条表示模拟值，而红色线条表示理论值（其中把概率理论值 的最大值调整为模拟值的最大值以便于对比）
(A)图源程序
clear; clc; close all; figure; xlim([-1,9]); ylim([0,9]); axis equal; hold on; L=0.8; No=8;% level number for N=1:No; mN=No/2-[1+N]/2; for k=1:N; plot([mN+k]*[1,1],[0,0.7]+No-N,'k','linewidth',2); text(-0.9,No-N+0.3,num2str(N),'fontsize',14,... 'fontname','times new roman'); end end arg=linspace(0,pi*2,200); col=[0.4,0.4,0.4]; fill(No/2+0.2*cos(arg),No+0.5+0.2*sin(arg),col,... 'Edgecolor',col); plot([No/2-1,No/2-1,-0.5,-0.5],... [No-1+0.7,No-1,0.7,0],'k'); plot([No/2+1,No/2+1,No+0.5,No+0.5],... [No-1+0.7,No-1,0.7,0],'k'); plot([-1,No+1],[0,0],'k','linewidth',2); for k=1:No+1; text(k-1.1,0.4,num2str(k),'fontsize',14,... 'fontname','times new roman'); plot(k-1+0.24*cos(arg),0.4+0.24*sin(arg),'k'); end set(gcf,'color','w'); axis off;
用Monte Carlo方法计算积分(随机投点法)
模拟赌博过程思路：在每一次模拟中,随机产生一个数,如果该数小于
clc; clear; close all; a=10; b=3; p=0.55; S=0; N=10000; m=6; %甲的赌本 %乙的赌本 %甲赢的概率 % 计数设置为0 % 模拟次数 %设定随机数状态值（1 2 3 4 5 6 ），改变这个值可以进行不同的实验
Monte Carlo仿真原理
Monte Carlo方法的的基本思想是首先建立一个概率模型，使所 求问题的解正好是该模型参数或其他有关特征量，然后通过模 拟 ———— 统计试验，即多次随机抽样试验（确定m和n）,统计 出某事件发生的百分比，只要试验次数很大，该百分比就近拟于 事件发生的概率。这实际上就是事件发生概率的统计定义。利用 建立的概率模型，求出要估计的参数。蒙特卡洛属于试验数学的 一个分支。