随机过程matlab程序.docx
- 文档编号:10313235
- 上传时间:2023-02-10
- 格式:DOCX
- 页数:26
- 大小:111.88KB
随机过程matlab程序.docx
《随机过程matlab程序.docx》由会员分享,可在线阅读,更多相关《随机过程matlab程序.docx(26页珍藏版)》请在冰豆网上搜索。
随机过程matlab程序
基本操作
-5/(4.8+5.32)^2
area=pi*2.5^2
x1=1+1/2+1/3+1/4+1/5+1/6
exp(acos(0.3))
a=[123;456;789]
a=[1:
3,4:
6,7:
9]
a1=[6:
-1:
1]
a=eye(4)a1=eye(2,3)b=zeros(2,10)c=ones(2,10)c1=8*ones(3,5)
d=zeros(3,2,2);
r1=rand(2,3)
r2=5-10*rand(2,3)
r4=2*randn(2,3)+3
arr1=[1.1-2.23.3-4.45.5]
arr1(3)arr1([14])arr1(1:
2:
5)
arr2=[123;-2-3-4;345]
arr2(1,:
)
arr2(:
1:
2:
3)
arr3=[12345678]
arr3(5:
end)arr3(end)
绘图
x=[0:
1:
10];
y=x.^2-10*x+15;
plot(x,y)
x=0:
pi/20:
2*pi
y1=sin(x);y2=cos(x);
plot(x,y1,'b-');
holdon;
plot(x,y2,‘k--’);
legend(‘sinx’,‘cosx’);
x=0:
pi/20:
2*pi;
y=sin(x);
figure
(1)
plot(x,y,'r-')
gridon
以二元函数图z=xexp(-x^2-y^2)为例讲解基本操作,首先需要利用meshgrid函数生成X-Y平面的网格数据,如下所示:
xa=-2:
0.2:
2;
ya=xa;
[x,y]=meshgrid(xa,ya);
z=x.*exp(-x.^2-y.^2);
mesh(x,y,z);
建立M文件
functionfenshu(grade)
ifgrade>95.0
disp('ThegradeisA.');
else
ifgrade>86.0
disp('ThegradeisB.');
else
ifgrade>76.0
disp('ThegradeisC.');
else
ifgrade>66.0
disp('ThegradeisD.');
else
disp('ThegradeisF.');
end
end
end
end
end
functiony=func(x)
ifabs(x)<1
y=sqrt(1-x^2);
elsey=x^2-1;
end
functionsumm(n)
i=1;
sum=0;
while(i<=n)
sum=sum+i;
i=i+1;
end
str=['¼ÆËã½á¹ûΪ£º',num2str(sum)];
disp(str)
end
求极限
symsx
limit((1+x)^(1/x),x,0,'right')
求导数
symsx;
f=(sin(x)/x);
diff(f)
diff(log(sin(x)))
求积分
symsx;
int(x^2*log(x))
symsx;
int(abs(x-1),0,2)
常微分方程求解
dsolve('Dy+2*x*y=x*exp(-x^2)','x')
计算偏导数
x/(x^2+y^2+z^2)^(1/2)
diff((x^2+y^2+z^2)^(1/2),x,2)
重积分
int(int(x*y,y,2*x,x^2+1),x,0,1)
级数
symsn;
symsum(1/2^n,1,inf)
Taylor展开式
求y=exp(x)在x=0处的5阶Taylor展开式
taylor(exp(x),0,6)
矩阵求逆
A=[0-6-1;62-16;-520-10]
det(A)
inv(A)
特征值、特征向量和特征多项式
A=[0-6-1;62-16;-520-10];
lambda=eig(A)
[v,d]=eig(A)
poly(A)
多项式的根与计算
p=[10-2-5];
r=roots(p)
p2=poly(r)
y1=polyval(p,4)
例子:
x=[-3:
3]'
y=[3.03,3.90,4.35,4.50,4.40,4.02,3.26]';
A=[2*x,2*y,ones(size(x))];
B=x.^2+y.^2;
c=inv(A'*A)*A'*B;
r=sqrt(c(3)+c
(1)^2+c
(2)^2)
例子
ezplot('-2/3*exp(-t)+5/3*exp(2*t)','-2/3*exp(-t)+2/3*exp(2*t)',[0,1])
gridon;axis([0,12,0,5])
密度函数和概率分布
设x~b(20,0.1),
binopdf(2,20,0.1)
分布函数
设x~N(1100,502),y~N(1150,802),则有
normcdf(1000,1100,50)=0.0228,1-0.0228=0.9772
normcdf(1000,1150,80)=0.0304,1-0.0304=0.9696
统计量数字特征
x=[29.827.628.3]
mean(x)
max(x)
min(x)
std(x)
symspk;
Ex=symsum(k*p*(1-p)^(k-1),k,1,inf)
symsxy;
f=x+y;
Ex=int(int(x*y*f,y,0,1),0,1)
参数估计
例:
对某型号的20辆汽车记录其5L汽油的行驶里程(公里),
观测数据如下:
29.8 27.6 28.3 27.9 30.1 28.7 29.9 28.0 27.9 28.7
28.4 27.2 29.5 28.5 28.0 30.0 29.1 29.8 29.6 26.9
设行驶里程服从正态分布,试用最大似然估计法求总体的均值和方差。
x1=[29.827.628.327.930.128.729.928.027.928.7];
x2=[28.427.229.528.528.030.029.129.829.626.9];
x=[x1x2]';
p=mle('norm',x);
muhatmle=p
(1),
sigma2hatmle=p
(2)^2
[m,s,mci,sci]=normfit(x,0.5)
假设检验
例:
下面列出的是某工厂随机选取的20只零部件的装配时间(分):
9.810.410.69.69.79.910.911.19.610.2
10.39.69.911.210.69.810.510.110.59.7
设装配时间总体服从正态分布,标准差为0.4,是否认定装配时间的均值在0.05的水平下不小于10。
解:
:
在正态总体的方差已知时MATLAB均值检验程序:
x1=[9.810.410.69.69.79.910.911.19.610.2];
x2=[10.39.69.911.210.69.810.510.110.59.7];
x=[x1x2]';m=10;sigma=0.4;a=0.05;[h,sig,muci]=ztest(x,m,sigma,a,1)
得到:
h=1,sig=0.873,muci=10.398Inf
%PPT例2一维正态密度与二维正态密度
symsxy;
s=1;t=2;
mu1=0;mu2=0;sigma1=sqrt((s^2));sigma2=sqrt((t^2));
x=-6:
0.1:
6;
f1=1/sqrt(2*pi*sigma1)*exp(-(x-mu1).^2/(2*sigma1^2));
f2=1/sqrt(2*pi*sigma2)*exp(-(x-mu2).^2/(2*sigma2^2));
plot(x,f1,'r-',x,f2,'k-.')
rho=(1+s*t)/(sigma1*sigma2);
f=1/(2*pi*sigma1*sigma2*sqrt(1-rho^2))*exp(-1/(2*(1-rho^2))*((x-mu1)^2/sigma1^2-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)^2/sigma2^2));
ezsurf(f)
%P34例3.1.1
p1=poisscdf(5,10)
p2=poisspdf(0,10)
[p1,p2]
%输出
p1=0.0671
p2=4.5400e-005
ans=0.06710.0000
%P40例3.2.1
p3=poisspdf(9,12)
%输出
p3=0.0874
%P40例3.2.2
p4=poisspdf(0,12)
%输出
p4=6.1442e-006
%P35-37(Th3.1.1)解微分方程
%输入:
symsp0p1p2;
S=dsolve('Dp0=-lamda*p0','Dp1=-lamda*p1+lamda*p0','Dp2=-lamda*p2+lamda*p1',
'p0(0)=1','p1(0)=0','p2(0)=0');
[S.p0,S.p1,S.p2]
%输出:
ans=
[exp(-lamda*t),exp(-lamda*t)*t*lamda,1/2*exp(-lamda*t)*t^2*lamda^2]
%P40泊松过程仿真
%simulate10times
clear;
m=10;lamda=1;x=[];
fori=1:
m
s=exprnd(lamda,'seed',1);
%seed是用来控制生成随机数的种子,使得生成随机数的个数是一样的.
x=[x,exprnd(lamda)];
t1=cumsum(x);
end
[x',t1']
%输出:
ans=
0.65090.6509
2.40613.0570
0.10023.1572
0.12293.2800
0.82334.1033
0.24634.3496
1.90746.2570
0.47836.7353
1.34478.0800
0.80828.8882
%输入:
N=[];
fort=0:
0.1:
(t1(m)+1)
ift (1) N=[N,0]; elseift (2) N=[N,1]; elseift N=[N,2]; elseift N=[N,3]; elseift N=[N,4]; elseift N=[N,5]; elseift N=[N,6]; elseift N=[N,7]; elseift N=[N,8]; elseift N=[N,9]; else N=[N,10]; end end plot(0: 0.1: (t1(m)+1),N,'r-') %输出: %simulate100times clear; m=100;lamda=1;x=[]; fori=1: m s=rand('seed'); x=[x,exprnd(lamda)]; t1=cumsum(x); end [x',t1'] N=[]; fort=0: 0.1: (t1(m)+1) ift (1) N=[N,0]; end fori=1: (m-1) ift>=t1(i)&t N=[N,i]; end end ift>t1(m) N=[N,m]; end end plot(0: 0.1: (t1(m)+1),N,'r-') %输出: %P48非齐次泊松过程仿真 %simulate10times clear; m=10;lamda=1;x=[]; fori=1: m s=rand('seed');%exprnd(lamda,'seed',1);setseeds x=[x,exprnd(lamda)]; t1=cumsum(x); end [x',t1'] N=[];T=[]; fort=0: 0.1: (t1(m)+1) T=[T,t.^3];%timeisadjusted,cumulativeintensityfunctionist^3. ift (1) N=[N,0]; end fori=1: (m-1) ift>=t1(i)&t N=[N,i]; end end ift>t1(m) N=[N,m]; end end plot(T,N,'r-') %output ans= 0.42200.4220 3.33233.7543 0.16353.9178 0.06833.9861 0.38754.3736 0.27744.6510 0.29694.9479 0.93595.8838 0.42246.3062 1.76508.0712 10timessimulation100timessimulation %P50复合泊松过程仿真 %simulate100times clear; niter=100;%iteratenumber lamda=1;%arrivingrate t=input('Inputatime: ','s') fori=1: niter rand('state',sum(clock)); x=exprnd(lamda);%intervaltime t1=x; whilet1 x=[x,exprnd(lamda)]; t1=sum(x);%arrivingtime end t1=cumsum(x); y=trnd(4,1,length(t1));%rand(1,length(t1));gamrnd(1,1/2,1,length(t1)));frnd(2,10,1,length(t1))); t2=cumsum(y); end [x',t1',y',t2'] X=[];m=length(t1); fort=0: 0.1: (t1(m)+1) ift (1) X=[X,0]; end fori=1: (m-1) ift>=t1(i)&t X=[X,t2(i)]; end end ift>t1(m) X=[X,t2(m)]; end end plot(0: 0.1: (t1(m)+1),X,'r-') 跳跃度服从[0,1]均匀分布情形跳跃度服从 分布情形 跳跃度服从t(10)分布情形 %%Simulatetheprobabilitythatsalesrevenuefallsinsomeinterval.(e.g.example3.3.6inteachingmaterial) clear; niter=1.0E4;%numberofiterations lamda=6;%arrivingrate(unit: minute) t=720;%12hours=720minutes above=repmat(0,1,niter);%setupstorage fori=1: niter rand('state',sum(clock)); x=exprnd(lamda);%intervaltime n=1; whilex x=x+exprnd(1/lamda);%arrivingtime ifx>=t n=n; else n=n+1; end end z=binornd(200,0.5,1,n);%generatensales y=sum(z); above(i)=sum(y>432000); end pro=mean(above) Output: pro=0.3192 %%Simulatethelosspro.ForaCompoundPoissonprocess clear; niter=1.0E3;%numberofiterations lamda=1;%arrivingrate t=input('Inputatime: ','s') below=repmat(0,1,niter);%setupstorage fori=1: niter rand('state',sum(clock)); x=exprnd(lamda);%intervaltime n=1; whilex x=x+exprnd(lamda);%arrivingtime ifx>=t n=n; else n=n+1; end end r=normrnd(0.05/253,0.23/sqrt(253),1,n);%generatenrandomjumps y=log(1.0E6)+cumsum(r); minX=min(y);%minmumreturnovernextnjumps below(i)=sum(minX end pro=mean(below) Output: t=50,pro=0.45 %P75(Example5.1.5)马氏链 chushivec0=[001000] P=[0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0] jueduivec1=chushivec0*P jueduivec2=chushivec0*(P^2) %计算1到n步后的分布 chushivec0=[001000]; P=[0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0]; n=10 t=1/6*ones([16]); jueduivec=repmat(t,[n1]); fork=1: n jueduiveck=chushivec0*(P^k); jueduivec(k,1: 6)=jueduiveck end %比较相邻的两行 n=70; jueduivecn=chushivec0*(P^n) n=71; jueduivecn=chushivec0*(P^n) %Replacethefirstdistribution,Comparingtwoneighbourabsolute-vectorsoncemore chushivec0=[1/61/61/61/61/61/6]; P=[0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0]; n=70; jueduivecn=chushivec0*(P^n) n=71; jueduivecn=chushivec0*(P^n) %赌博问题模拟(带吸收壁的随机游走: 结束1次游走所花的时间及终止状态) a=5;p=1/2; m=0; whilem<100 m=m+1; r=2*binornd(1,p)-1; ifr==-1 a=a-1; else a=a+1; end ifa==0|a==10 break; end end [ma] %赌博问题模拟(带吸收壁的随机游走: 结束N次游走所花的平均时间及终止状态分布规律) %p=q=1/2 p=1/2; m1=0;m2=0;N=1000; t1=0;t2=0; forn=1: 1: N m=0;a=5; whilea>0&a<10 m=m+1; r=2*binornd(1,p)-1; ifr==-1 a=a-1; else a=a+1; end end ifa==0 t1=t1+m;m1=m1+1; else t2=t2+m;m2=m2+1; end end fprintf('Theaveragetimesofarriving0and10respectivelyare%d,%d.\n',[t1/m1,t2/m2]); fprintf('Thefrequenciesofarriving0and10respectivelyare%d,%d.\n',[m1/N,m2/N]); %verify: fprintf('Theprobabilityofarriving0anditsapproximaterespectivelyare%d,%d.\n',[5/10,m1/N]); fprintf('Theexpectationofarriving0or10anditsapproximaterespectivelyare%d,%d.\n',[5*(10-5)/(2*p),(t1+t2)/N]); %p~=q p=1/4; m1=0;m2=0;N=1000; t1=zeros(1,N);t2=zeros(1,N); forn=1: 1: N m=0;a=5; whilea>0&a<15 m=m+1; r=2*binornd(1,p)-1; ifr==-1 a=a-1; else a=a+1; end end ifa==0 t1(1,n)=m;m1=m1+1; else t2(1,n)=m;m2=m2+1; end end fprintf('T
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 随机 过程 matlab 程序
![提示](https://static.bdocx.com/images/bang_tan.gif)