矩阵运算1）矩阵输入： A=13 -56 78;25 63 -235;78 25 563;1 0 -1或或A=13,-56,78;25,63,-235;78,25,563;1,0,-12）矩阵的运算A+B; A-B; c*A ; A*B; A ; A逆: inv(A);矩阵除法：AB等效于inv(A)*B， B/A等效于 B*inv(A)。方阵行列式值det(A);矩阵的秩rank(A)矩阵的迹trace(A)7/27/20241用矩阵求解线性方程组线性方程组表示为Ax=b，A不一定是方阵。命令格式：C=A b; rref(C)例：求解线性方程组 x+2y+3z+t=3 x+4y+5z+2t=2 2x+9y+8z+3t=7 3x+7y+7z+2t=12 命令：A=1 2 3 1;1 4 5 2;2 9 8 3;3 7 7 2;b=3;2;7;12;C=A b; rref(C)7/27/20242一、基本实验命令：1、矩阵输入及运算: 2、矩阵的乘方: An 3 、 求解平稳分布归结为求线性方程组表示为Ax=b，A不一定是方阵。命令格式：C=A b; rref(C)实验实验1 验证马氏链的遍历性验证马氏链的遍历性 及求平稳分布及求平稳分布7/27/20243题目、求解线性方程组x+2y+3z+t=3 x+4y+5z+2t=2 2x+9y+8z+3t=7 3x+7y+7z+2t=12命令：A=1 2 3 1;1 4 5 2;2 9 8 3;3 7 7 2;b=3;2;7;12;C=A b; rref(C)7/27/20244A=1 2 3 1;1 4 5 2;2 9 8 3;3 7 7 2;b=3;2;7;12;C=A b; rref(C)ans = 1.0000 0 0 -0.5000 5.1667 0 1.0000 0 0 0.6667 0 0 1.0000 0.5000 -1.1667 0 0 0 0 0即x-0.5t=5.1667y=0.6667z+0.5t=-1.16677/27/20245基本原理：基本原理：则此马氏链是遍历的，即是方程组j =0，1，2，s的唯一解平稳分布平稳分布二、提高实验7/27/20246其一步转移矩阵为试证此链具有遍历性，并求出平稳分布。题目题目1解解二、提高实验7/27/20247所以因此，该马氏链具有遍历性。解得所以马氏链的平稳分布为X123各状态的平均返回时间各状态的平均返回时间7/27/20248命令：P=1/3 2/3 0;1/3 0 2/3 ;0 1/3 2/3;PnB =1/3 2/3 0 1;1/3 0 2/3 1;0 1/3 2/3 1;D=1 0 0 0;0 1 0 0;0 0 1 0; A=B-D;or A=-1+1/3 1/3 0;2/3 -1 1/3 ;0 2/3 -1+2/3;1 1 1;b=0;0;0;1;C=A b; rref(C)7/27/20249其一步转移矩阵为试证此链具有遍历性，并求出平稳分布，平均返回时间。题目题目2 7/27/202410所以因此，该马氏链具有遍历性。解得7/27/202411市场占有率预测市场占有率预测设某地有设某地有1600户居民，某产品只有甲、乙、丙户居民，某产品只有甲、乙、丙3厂厂家在该地销售。经调查，家在该地销售。经调查，8月份买甲、乙、丙三厂月份买甲、乙、丙三厂的户数分别为的户数分别为480，320，800。9月份里，原买甲月份里，原买甲的有的有48户转买乙产品，有户转买乙产品，有96户转买丙产品；原买户转买丙产品；原买乙的有乙的有32户转买甲产品，有户转买甲产品，有64户转买丙产品；原户转买丙产品；原买丙的有买丙的有64户转买甲产品，有户转买甲产品，有32户转买乙产品。户转买乙产品。用状态用状态1、2、3分别表示甲、乙、丙三厂，试求分别表示甲、乙、丙三厂，试求（1）转移概率矩阵；）转移概率矩阵；（2）9月份市场占有率的分布；月份市场占有率的分布；（3）12月份市场占有率的分布；月份市场占有率的分布；（4）当顾客流如此长期稳定下去市场占有率的分布。）当顾客流如此长期稳定下去市场占有率的分布。（5）各状态的平均返回时间）各状态的平均返回时间实验题目实验题目三、创新实验7/27/202412解解（1） 由题意得频数转移矩阵为由题意得频数转移矩阵为再用频数估计概率，得转移概率矩阵为（2）以）以1600除以除以N中各行元素之和，中各行元素之和，得初始概率分布（即初始市场占有率）得初始概率分布（即初始市场占有率）7/27/202413所以9月份市场占有率分布为（3）12月份市场占有率分布为7/27/202414（4）由于该链不可约、非周期、状态有限正常返的，所以是遍历的。解方程组即得当顾客流如此长期稳定下去是市场占有率的分布为7/27/202415命令： p=336/480 48/480 96/480;32/320 224/320 64/320;64/800 32/800 704/800;p0=0.3 0.2 0.5;p1=p0*pp4= p0*p4B =336/480 48/480 96/480 1;32/320 224/320 64/320 1;64/800 32/800 704/800 1;D=1 0 0 0;0 1 0 0;0 0 1 0; A=B-D;b=0;0;0;1;C=A b; rref(C)7/27/202416下面数据为1915到1973年的最大年径流量（按行从左到右）：15600 8960 10400 10600 10800 9880 9850 10900 8810 9960 12200 7510 8640 6380 6810 8820 14400 7440 7240 6430 11000 7340 9260 5290 9130 7480 6980 9650 7260 8750 9900 7310 9040 7310 8850 7840 10700 6190 9610 7580 9990 6150 8250 6030 8980 6180 9630 9490 2340 11100 5090 10900 6490 12600 6640 7430 6760 10000 93007/27/202417实验实验2 时间序列的线形模型的识别时间序列的线形模型的识别一、基本实验一、基本实验1、平稳平稳时间序列的标准化时间序列的标准化z=15600 8960 10400 10600 10800 9880 9850 10900 8810 9960 12200 7510 +.8640 6380 6810 8820 14400 7440 7240 6430 11000 7340 9260 5290 +.9130 7480 6980 9650 7260 8750 9900 7310 9040 7310 8850 7840 +.10700 6190 9610 7580 9990 6150 8250 6030 8980 6180 9630 9490 +.2340 11100 5090 10900 6490 12600 6640 7430 6760 10000 9300;命令：（1） %计算z的平均值 m=mean(z) m= 8.6686e+003 （2） %z的标准化的标准化 for k=1:59 w(k)=z(k)-m; end w7/27/2024187/27/202419二、提高实验基本原理：基本原理： 模型 AR(p) MA(q) ARMA(p,q) 拖尾 截尾在 k=q处 拖尾 截尾在 k=p处 拖尾 拖尾 7/27/2024201、自相关函数及图形命令： %求X的自相关函数（1）ACF=autocorr(X,13)ACF = 1.0000 -0.2305 0.2928 -0.1628 0.2831 -0.0090 0.2186 0.0774 -0.0042 0.0521 0.0247 0.0856 -0.0738 0.0313 %画出自相关函数图形 （2） autocorr(X,13) 7/27/202421画出自相关函数p的图像画出自相关函数p的图像7/27/2024222、偏相关函数及图形命令：（1） %求X的偏相关函数值PartialACF=parcorr(X,13)PartialACF = 1.0000 -0.2310 0.2422 -0.0517 0.1910 0.1875 0.2024 0.2792 0.0079 0.0756 0.0690 0.1197 0.0085 -0.0014（2）%画出偏相关函数的图形 parcorr(X,13)7/27/202423画出偏相关函数7/27/202424%同时画出自相关函数图形和偏相关函数的图形命令：subplot(121);autocorr(X,13) subplot(122);parcorr(X,13)7/27/2024257/27/202426三、创新实验某种股票的价格数据如下：某种股票的价格数据如下：z=116 90 104 106 108 99 98 109 88 99 122 75 86 64 68 88 114 74 72 64 110 73 92 52 91 74 69 96 72 87 99 73 90 73 88 78 107 61 96 75 99 61 82 60 89 ;判别该数据符合哪类线性模型？判别该数据符合哪类线性模型？命令：（1）%验证股票的价格数据是平稳时间序列股票的价格数据是平稳时间序列 plot(z)（2） %计算z的平均值 m=mean(z)（3） %z的标准化的标准化 for k=1:45 X(k)=z(k)-m; end X7/27/2024277/27/202428m = 86.4667X= Columns 1 through 20 29.5333 3.5333 17.5333 19.5333 21.5333 12.5333 11.5333 22.5333 1.5333 12.5333 35.5333 -11.4667 -0.4667 -22.4667 -18.4667 1.5333 27.5333 -12.4667 -14.4667 -22.4667 Columns 21 through 40 23.5333 -13.4667 5.5333 -34.4667 4.5333 -12.4667 -17.4667 9.5333 -14.4667 0.5333 12.5333 -13.4667 3.5333 -13.4667 1.5333 -8.4667 20.5333 -25.4667 9.5333 -11.4667 Columns 41 through 45 12.5333 -25.4667 -4.4667 -26.4667 2.53337/27/2024291、自相关函数及图形命令： %求w的自相关函数（1）ACF=autocorr(X,10)ACF = 1.0000 -0.5684 0.3709 -0.3144 0.2685 -0.2082 0.1240 0.0658 -0.0097 -0.0504 -0.0292 %画出自相关函数图形 （2） autocorr(X,10) 7/27/202430画出自相关函数p的图像画出自相关函数p的图像7/27/2024312、偏相关函数及图形命令：（1） %求w的偏相关函数值PartialACF=parcorr(X,10)PartialACF = 1.0000 -0.5716 0.0619 -0.1197 0.0681 0.0186 -0.0975 0.3536 0.2211 -0.2127 -0.6025（2）%画出偏相关函数的图形 parcorr(X,10)7/27/202432画出偏相关函数7/27/202433%同时画出自相关函数图形和偏相关函数的图形命令：subplot(121);autocorr(X,10) subplot(122);parcorr(X,10)自相关函数截尾；偏相关函数拖尾。是MA（q）模型7/27/2024347/27/202435实验实验2 时间序列的线形模型的识别时间序列的线形模型的识别1、平稳平稳时间序列的标准化时间序列的标准化z=15600 8960 10400 10600 10800 9880 9850 10900 8810 9960 12200 7510 8640 6380 6810 8820 14400 7440 7240 6430 11000 7340 +. 9260 5290 9130 7480 6980 9650 7260 8750 9900 7310 9040 7310 +. 8850 7840 10700 6190 9610 7580 9990 6150 8250 6030 8980 +.6180 9630 9490 2340 11100 5090 10900 6490 12600 6640 +. 7430 6760 10000 9300;命令：m=mean(z) %计算z的平均值 for k=1:59 X(k)=z(k)-m; %z的标准化的标准化 end X7/27/2024362、求自相关函数p% 求自相关函数prc=0;for k=1:59rc=rc+1/59*X(k)*X(k);end% rc 为r的初始值for k=1:15x=0;for j=1:59-kx=x+1/59*X(j)*X(j+k);endr(k)=x;p(k)=r(k)/rc;end7/27/202437画出自相关函数p的图像%画出自相关函数p的图像x=1:15;plot(x,p,x,0,-k,x,p,*); xlabel(k的取值/k)ylabel(自相关函数值 /p)画出自相关函数p的图像7/27/2024383、求偏相关函数m(k,k)%求偏相关函数m(k,k)m(1,1)=p(1);for k=1:14x=0;y=0;for j=1:kx=x+p(k+1-j)*m(k,j);y=y+p(j)*m(k,j);endm(k+1,k+1)=(p(k+1)-x)/(1-y);for j=1:km(k+1,j)=m(k,j)-m(k+1,k+1)*m(k,k-j+1);endend(diag(m)计算出的偏相关函数m(k,k)如下7/27/202439画出偏相关函数m(k,k)(取k=15)图像接着在MATLAB命令窗口中输入下列程序：% 画出偏相关函数m(k,k)的图像x=1:15;plot(x,(diag(m),x,0,-k,x,(diag(m), *) xlabel(k的取值)ylabel (偏相关函数值 /m)7/27/202440偏相关函数m(k,k)(取k=15)图像7/27/202441实验实验3 平稳时间序列线性模型的计算及预测平稳时间序列线性模型的计算及预测一、基本实验一、基本实验1、命令：a=lpc(z,p)Z是时间序列的数据，p表示AR(p)的阶数2、命令：estx=filter(0-a(2:end),1,z)%表示估计时间序列3、命令：e=z-estx%表示估计时间序列的预测误差4、命令:acs,lags=zcorr(e,coeff)%表示时间序列预测误差的自相关函数5、命令：plot(1:59,z,1:59,estx,-.)%表示时间序列原始数据与预测值的叠加图6、plot（lags,acs)%时间序列预测误差的自相关函数图7/27/202442例：例： z=15600 8960 10400 10600 10800 9880 9850 10900 8810 9960 12200 7510 +.8640 6380 6810 8820 14400 7440 7240 6430 11000 7340 9260 5290 +.9130 7480 6980 9650 7260 8750 9900 7310 9040 7310 8850 7840 +.10700 6190 9610 7580 9990 6150 8250 6030 8980 6180 9630 9490 +.2340 11100 5090 10900 6490 12600 6640 7430 6760 10000 9300; a=lpc(z,2)a = 1.0000 -0.4104 -0.5406 estx=filter(0-a(2:end),1,z)； estx = 1.0e+004 * Columns 1 through 16 0.6401 1.2110 0.9111 0.9972 1.0162 0.9893 0.9383 0.9798 0.9508 0.8850 1.0390 0.9677 0.7605 0.7289 0.6243 0.7301 Columns 17 through 32 1.0677 1.0837 0.6993 0.6552 0.7990 0.8958 0.7768 0.7177 0.6606 0.8005 0.6908 0.7733 0.8196 0.7515 0.8793 0.8351 Columns 33 through 48 0.7661 0.7887 0.7583 0.8001 0.8629 0.8324 0.7290 0.8305 0.8197 0.7924 0.6710 0.6934 0.6945 0.7390 0.7292 0.9100 Columns 49 through 59 0.6090 0.5820 0.8089 0.7224 0.8556 0.8679 0.9536 0.6638 0.6790 0.7758 0.9222 e=z-estx7/27/202443 plot(1:59,z,1:59,estx,-.);7/27/202444 acs,lags= xcorr(e,coeff); plot(lags,acs);7/27/202445plot(1:59,z,1:59,estx,-.); title(原始信号); xlabel(采样点);ylabel(幅度);grid; legend(原始信号,估计值)7/27/202446plot(lags,acs);title(预测误差的自相关函数); xlabel(延迟);ylabel(归一化值);grid;7/27/202447左右画出两张图subplot(121) plot(1:59,z,1:59,estx,-.); title(原始信号);xlabel(采样点);ylabel(幅度);grid;legend(原始信号,估计值)subplot(122)plot(lags,acs);title(预测误差的自相关函数);xlabel(延迟);ylabel(归一化值);grid;7/27/2024487/27/202449上下画两张图 subplot(211) plot(1:59,z,1:59,estx,-.); title(原始信号);%图上面的注释xlabel(采样点);ylabel(幅度);grid;legend(原始信号,估计值)subplot(212)plot(lags,acs);title(预测误差的自相关函数);xlabel(延迟);ylabel(归一化值);grid;7/27/2024507/27/202451二、提高实验基本原理：基本原理：1、确定、确定 AR（p）模型）模型 参数参数2、确定、确定MA（q） 模型参数模型参数3、确定、确定ARMA（p,q） 模型参数模型参数7/27/202452图形的处理方法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);7/27/202453三、创新实验已知某一地区1990-2009年的肿瘤引起的死亡率如下：利用数据建立AR（3）模型，并分析误差死亡率1011.3 99.19.4 9.18.78.6 9.1 9.210.17.36.6 5.75.9 5.65.55.6 5.7 6.47/27/202454z=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.6,5.6,5.7,6.4; % 求求AR(p)模型的参数 a=lpc(z,3)a= %求估计的时间序列值estx=filter(0-a(2:end),1,z) estx =%画出原数据与估计的时间序列值的叠加图plot(1:20,z,1:20,estx,-.); 7/27/202455%计算误差计算误差e=z-estxe=%画出误差自相关函数图acs,lags= xcorr(e,coeff); plot(lags,acs);%画出带注释的左右图形 subplot(121) plot(1:20,z,1:20,estx,-.); title(原始信号);xlabel(采样点);ylabel(幅度);grid;legend(原始信号,估计值) subplot(121) plot(lags,acs);title(预测误差的自相关函数);xlabel(延迟);ylabel(归一化值);grid;7/27/202456%画出带注释的上下图形 subplot(211) plot(1:20,z,1:20,estx,-.); title(原始信号);xlabel(采样点);ylabel(幅度);grid;legend(原始信号,估计值) subplot(212) plot(lags,acs);title(预测误差的自相关函数);xlabel(延迟);ylabel(归一化值);grid;7/27/202457