随机过程matlab程序Word文件下载.docx
《随机过程matlab程序Word文件下载.docx》由会员分享，可在线阅读，更多相关《随机过程matlab程序Word文件下载.docx（26页珍藏版）》请在冰点文库上搜索。
plot（x,y2,‘k--’）;
legend（‘sinx’,‘cosx’）;
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
ThegradeisB.'
else
76.0
ThegradeisC.'
66.0
ThegradeisD.'
ThegradeisF.'
end
end
functiony=func（x）
ifabs（x）<
1
y=sqrt（1-x^2）;
elsey=x^2-1;
functionsumm（n）
i=1;
sum=0;
while（i<
=n）
sum=sum+i;
i=i+1;
str=['
¼
Æ
Ë
ã
½
á
¹
û
Î
ª
£
º
'
num2str（sum）];
disp（str）
求极限
symsx
limit（（1+x）^（1/x）,x,0,'
right'
求导数
symsx;
f=（sin（x）/x）;
diff（f）
diff（log（sin（x）））
求积分
int（x^2*log（x））
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）
特征值、特征向量和特征多项式
-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];
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一维正态密度与二维正态密度
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,'
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;
lamda=1;
x=[];
fori=1:
m
s=exprnd（lamda,'
seed'
1）;
%seed是用来控制生成随机数的种子,使得生成随机数的个数是一样的.
x=[x,exprnd（lamda）];
t1=cumsum（x）;
[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:
（t1（m）+1）
ift<
t1
（1）
N=[N,0];
elseift<
t1
（2）
N=[N,1];
t1（3）
N=[N,2];
t1（4）
N=[N,3];
t1（5）
N=[N,4];
t1（6）
N=[N,5];
t1（7）
N=[N,6];
t1（8）
N=[N,7];
t1（9）
N=[N,8];
t1（10）
N=[N,9];
N=[N,10];
plot（0:
（t1（m）+1）,N,'
%simulate100times
m=100;
s=rand（'
（m-1）
ift>
=t1（i）&
t<
t1（i+1）
N=[N,i];
ift>
t1（m）
N=[N,m];
%P48非齐次泊松过程仿真
s=rand（'
%exprnd（lamda,'
setseeds
T=[];
T=[T,t.^3];
%timeisadjusted,cumulativeintensityfunctionist^3.
plot（T,N,'
%output
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复合泊松过程仿真
niter=100;
%iteratenumber
lamda=1;
%arrivingrate
t=input（'
Inputatime:
s'
niter
rand（'
state'
sum（clock））;
x=exprnd（lamda）;
%intervaltime
t1=x;
whilet1<
t
x=[x,exprnd（lamda）];
t1=sum（x）;
%arrivingtime
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）;
y'
t2'
X=[];
m=length（t1）;
ift<
X=[X,0];
fori=1:
X=[X,t2（i）];
X=[X,t2（m）];
（t1（m）+1）,X,'
跳跃度服从[0,1]均匀分布情形跳跃度服从
分布情形
跳跃度服从t（10）分布情形
%%Simulatetheprobabilitythatsalesrevenuefallsinsomeinterval.（e.g.example3.3.6inteachingmaterial）
niter=1.0E4;
%numberofiterations
lamda=6;
%arrivingrate（unit:
minute）
t=720;
%12hours=720minutes
above=repmat（0,1,niter）;
%setupstorage
n=1;
whilex<
t
x=x+exprnd（1/lamda）;
%arrivingtime
ifx>
=t
n=n;
else
n=n+1;
end
z=binornd（200,0.5,1,n）;
%generatensales
y=sum（z）;
above（i）=sum（y>
432000）;
pro=mean（above）
Output:
pro=0.3192
%%Simulatethelosspro.ForaCompoundPoissonprocess
niter=1.0E3;
below=repmat（0,1,niter）;
%setupstorage
x=x+exprnd（lamda）;
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<
log（950000））;
pro=mean（below）
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];
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
%比较相邻的两行
n=70;
jueduivecn=chushivec0*（P^n）
n=71;
%Replacethefirstdistribution,Comparingtwoneighbourabsolute-vectorsoncemore
chushivec0=[1/61/61/61/61/61/6];
%赌博问题模拟（带吸收壁的随机游走：
结束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;
ifa==0|a==10
break;
[ma]
结束N次游走所花的平均时间及终止状态分布规律）
%p=q=1/2
m1=0;
m2=0;
N=1000;
t1=0;
t2=0;
forn=1:
N
m=0;
a=5;
whilea>
0&
a<
10
ifa==0
t1=t1+m;
m1=m1+1;
t2=t2+m;
m2=m2+1;
fprintf（'
Theaveragetimesofarriving0and10respectivelyare%d,%d.\n'
[t1/m1,t2/m2]）;
Thefrequenciesofarriving0and10respectivelyare%d,%d.\n'
[m1/N,m2/N]）;
%verify:
Theprobabilityofarriving0anditsapproximaterespectivelyare%d,%d.\n'
[5/10,m1/N]）;
Theexpectationofarriving0or10anditsapproximaterespectivelyare%d,%d.\n'
[5*（10-5）/（2*p）,（t1+t2）/N]）;
%p~=q
p=1/4;
t1=zeros（1,N）;
t2=zeros（1,N）;
15
t1（1,n）=m;
t2（1,n）=m;
T
升级会员 