随机过程实验课件三次

2023/5/221一、基本实验命令：1、矩阵输入及运算:2、矩阵的乘方:A^n3、求解平稳分布归结为求线性方程组表示为Ax=b，A不一定是方阵。命令格式：C=[Ab];rref(C)实验1验证马氏链的遍历性

及平稳分布2023/5/222【理论总结】遍历性平稳分布判别定理平稳分布计算2023/5/223例：求解线性方程组x+2y+3z+t=3x+4y+5z+2t=22x+9y+8z+3t=73x+7y+7z+2t=12命令：A=[1231;1452;2983;3772];b=[3;2;7;12];C=[Ab];rref(C)

复习：解线性方程组的实验命令2023/5/224A=[1231;1452;2983;3772];b=[3;2;7;12];C=[Ab];rref(C)ans=1.000000-0.50005.166701.0000000.6667001.00000.5000-1.166700000即x-0.5t=5.1667y=0.6667z+0.5t=-1.16672023/5/225基本原理：则此马氏链是遍历的，即是方程组j=0，1，2，…，s的唯一解平稳分布二、提高实验2023/5/226其一步转移矩阵为试证此链具有遍历性，并求出平稳分布。题目1解二、提高实验2023/5/227所以因此，该马氏链具有遍历性。解得所以马氏链的平稳分布为X123各状态的平均返回时间2023/5/228命令：P=[1/32/30;1/302/3;01/32/3];P^nB=[1/32/301;1/302/31;01/32/31];D=[1000;0100;0010];A=B’-D’;orA=[-1+1/31/30;2/3-11/3;02/3-1+2/3;111];b=[0;0;0;1];C=[Ab];rref(C)2023/5/229其一步转移矩阵为试证此链具有遍历性，并求出平稳分布，平均返回时间。首页题目2

2023/5/2210所以因此，该马氏链具有遍历性。解得2023/5/2211市场占有率预测设某地有1600户居民，某产品只有甲、乙、丙3厂家在该地销售。经调查，8月份买甲、乙、丙三厂的户数分别为480，320，800。9月份里，原买甲的有48户转买乙产品，有96户转买丙产品；原买乙的有32户转买甲产品，有64户转买丙产品；原买丙的有64户转买甲产品，有32户转买乙产品。用状态1、2、3分别表示甲、乙、丙三厂，试求（1）转移概率矩阵；（2）9月份市场占有率的分布；（3）12月份市场占有率的分布；（4）当顾客流如此长期稳定下去市场占有率的分布。（5）各状态的平均返回时间实验题目三、创新实验2023/5/2212解（1）由题意得频数转移矩阵为再用频数估计概率，得转移概率矩阵为（2）以1600除以N中各行元素之和，得初始概率分布（即初始市场占有率）2023/5/2213所以9月份市场占有率分布为（3）12月份市场占有率分布为2023/5/2214（4）由于该链不可约、非周期、状态有限正常返的，所以是遍历的。解方程组即得当顾客流如此长期稳定下去是市场占有率的分布为返回首页2023/5/2215命令：p=[];p0=[];p1=p0*pp4=p0*p^4B=[];D=[];A=B'-D';b=[0;0;0;1];C=[Ab];rref(C)2023/5/22162023/5/22172023/5/2218

B=[3/51/51/101/101;1/57/101/1001;001/21/21;1/3002/31];D=[10000;01000;00100;00010];A=B’-D’;b=[0;0;0;0;1];C=[Ab];rref(C)ans=1.00000000.357101.0000000.2381001.000000.11900001.00000.285700000

地层剖面岩性分布计算的Matlab程序2023/5/2219下面数据为1915到1973年的最大年径流量（按行从左到右）：1560089601040010600108009880985010900881099601220075108640638068108820144007440724064301100073409260529091307480698096507260875099007310904073108850784010700619096107580999061508250603089806180963094902340111005090109006490126006640743067601000093002023/5/2220实验2时间序列的线形模型的识别一、基本实验1、平稳时间序列的标准化z=[156008960104001060010800988098501090088109960122007510+...86406380681088201440074407240643011000734092605290+...913074806980965072608750990073109040731088507840+...1070061909610758099906150825060308980618096309490+...234011100509010900649012600664074306760100009300];命令：（1）%计算z的平均值

m=mean(z)m=8.6686e+003

（2）%z的标准化

fork=1:59X(k)=z(k)-m;endX2023/5/22212023/5/2222二、提高实验基本原理：

模型AR(p)MA(q)ARMA(p,q)

拖尾

截尾在k=q处

拖尾

截尾在k=p处

拖尾

拖尾2023/5/2223第一步：自相关函数及图形命令：%求X的自相关函数（1）[ACF]=autocorr(X,13)ACF=1.0000-0.23050.2928-0.16280.2831-0.00900.21860.0774-0.00420.05210.02470.0856-0.07380.0313%画出自相关函数图形（2）autocorr(X,13)

2023/5/2224画出自相关函数p的图像画出自相关函数p的图像2023/5/2225第2步：偏相关函数及图形

命令：（1）

%求X的偏相关函数值[PartialACF]=parcorr(X,13)PartialACF=1.0000-0.23100.2422-0.05170.19100.18750.20240.27920.00790.07560.06900.11970.0085-0.0014（2）%画出偏相关函数的图形>>parcorr(X,13)2023/5/2226画出偏相关函数2023/5/2227%同时画出自相关函数图形和偏相关函数的图形命令：subplot(121);autocorr(X,13)subplot(122);parcorr(X,13)2023/5/22282023/5/2229三、创新实验

题目：某种股票的价格数据如下：z=[116901041061089998109889912275866468881147472641107392529174699672879973907388781076196759961826089];判别该数据符合哪类线性模型？命令：（1）%验证股票的价格数据是平稳时间序列

plot(z)（2）%计算z的平均值

m=mean(z)（3）%z的标准化

fork=1:45X(k)=z(k)-m;endX2023/5/22302023/5/2231m=86.4667X=Columns1through2029.53333.533317.533319.533321.533312.533311.533322.53331.533312.533335.5333-11.4667-0.4667-22.4667-18.46671.533327.5333-12.4667-14.4667-22.4667Columns21through4023.5333-13.46675.5333-34.46674.5333-12.4667-17.46679.5333-14.46670.533312.5333-13.46673.5333-13.46671.5333-8.466720.5333-25.46679.5333-11.4667Columns41through4512.5333-25.4667-4.4667-26.46672.53332023/5/22321、自相关函数及图形命令：%求w的自相关函数（1）[ACF]=autocorr(X,10)ACF=1.0000-0.56840.3709-0.31440.2685-0.20820.12400.0658-0.0097-0.0504-0.0292%画出自相关函数图形（2）autocorr(X,10)

2023/5/2233画出自相关函数p的图像画出自相关函数p的图像2023/5/22342、偏相关函数及图形

命令：（1）

%求w的偏相关函数值[PartialACF]=parcorr(X,10)PartialACF=

1.0000-0.57160.0619-0.11970.06810.0186-0.09750.35360.2211-0.2127-0.6025（2）%画出偏相关函数的图形>>parcorr(X,10)2023/5/2235画出偏相关函数2023/5/2236%同时画出自相关函数图形和偏相关函数的图形命令：subplot(121);autocorr(X,10)subplot(122);parcorr(X,10)自相关函数截尾；偏相关函数拖尾。是MA（q）模型2023/5/22372023/5/2238

实验3平稳时间序列线性模型的计算及预测

一、基本实验1、命令：a=lpc(x,p)%求AR(p)模型的参数Z是时间序列的数据，p表示AR(p)的阶数命令[b,a]=prony(x,q,p)%求ARMA(p,q)模型的参数2、命令：estx=filter([0-a(2:end),1,x)%表示估计时间序列3、命令：e=z-estx%表示估计时间序列的预测误差4、命令:[acs,lags]=zcorr(e,’coeff’)%表示时间序列预测误差的自相关函数5、命令：plot(1:59,x,1:59,estx,‘-.’)%表示时间序列原始数据与预测值的叠加图6、plot（lags,acs)%时间序列预测误差的自相关函数图2023/5/2239

例：

z=[156008960104001060010800988098501090088109960122007510+...

86406380681088201440074407240643011000734092605290+...

913074806980965072608750990073109040731088507840+...

1070061909610758099906150825060308980618096309490+...

234011100509010900649012600664074306760100009300];

>>a=lpc(z,2)

a=

1.0000-0.4104-0.5406

>>estx=filter([0-a(2:end)],1,z)；estx=1.0e+004*Columns1through160.64011.21100.91110.99721.01620.98930.93830.97980.95080.88501.03900.96770.76050.72890.62430.7301Columns17through321.06771.08370.69930.65520.79900.89580.77680.71770.66060.80050.69080.77330.81960.75150.87930.8351Columns33through480.76610.78870.75830.80010.86290.83240.72900.83050.81970.79240.67100.69340.69450.73900.72920.9100Columns49through590.60900.58200.80890.72240.85560.86790.95360.66380.67900.77580.9222>>e=z-estx

2023/5/2240

plot(1:59,z,1:59,estx,'-.');2023/5/2241

[acs,lags]=xcorr(e,'coeff');plot(lags,acs);

2023/5/2242plot(1:59,z,1:59,estx,'-.');

title('原始信号');

xlabel('采样点');ylabel('幅度');grid;

legend('原始信号','估计值')2023/5/2243plot(lags,acs);

title('预测误差的自相关函数');

xlabel('延迟');ylabel('归一化值');grid;2023/5/2244左右画出两张图subplot(121)plot(1:59,z,1:59,estx,'-.');title('原始信号');xlabel('采样点');ylabel('幅度');grid;legend('原始信号','估计值')subplot(122)plot(lags,acs);title('预测误差的自相关函数');xlabel('延迟');ylabel('归一化值');grid;2023/5/22452023/5/2246上下画两张图subplot(211)plot(1:59,z,1:59,estx,'-.');title(‘原始信号’);%图上面的注释xlabel('采样点');ylabel('幅度');grid;legend('原始信号','估计值')subplot(212)plot(lags,acs);title('预测误差的自相关函数');xlabel('延迟');ylabel('归一化值');grid;2023/5/22472023/5/2248二、提高实验基本原理：1、确定AR（p）模型参数2、确定MA（q）模型参数3、确定ARMA（p,q）模型参数2023/5/2249图形的处理方法plot(1:59,z,1:59,estx,'-.');

title('原始信号');

xlabel('采样点');ylabel('幅度');grid;

legend('原始信号','估计值')plot(lags,acs);

title('预测误差的自相关函数');

xlabel('延迟');ylabel('归一化值');grid;subplot(121)plot(1:59,z,1:59,estx,'-.');subplot(122)plot(lags,acs);subplot(211)plot(1:59,z,1:59,estx,'-.');subplot(212)plot(lags,acs);2023/5/2250三、创新实验已知某一地区1990-2009年的肿瘤引起的死亡率如下：利用数据建立AR（3）模型，并分析误差死亡率‰1011.399.19.49.18.78.69.19.210.17.36.65.75.95.65.55.65.76.42023/5/2251

z=[10,11.3,9,9.1,9.4,9.1,8.7,8.6,9.1,9.2,10.1,7.3,6.2,7.3,6.2,5.7,5