时间序列分析 2.docx

1、时间序列分析 2评分20142015学年度第 3 学期一、 实验目的在生物医学领域，时间序列分析被广泛用于生理信号分析（如心电信号、脑电信号）流行病学研究（如传染病发病情况的预测）医学图像特征提取、（如fMEI图像序列中功能活动的激活点提取）等方面。本次实验要求熟悉时间系列分析的移动平均法、滑动平均法、自相关函数分析法，了解脑电信号分析的matlab示例。二、实验环境 1、硬件配置：处理器：Intel（R） Pentium(R) CPU 1.80GHz 3.1880GHz 安装内存（RAM）：4.00GB 系统类型：64位操作系统 2、软件环境：MATLAB R2010b软件三、实验内容(包括

2、本实验要完成的实验问题及需要的相关知识简单概述)1. 实验原理趋势分析方法（1）移动平均法。 设某一时间序列为， 则时刻的预测值为 （2）滑动平均法 式中； 为t点的滑动平均值；为单侧平滑时距 （3）指数平滑法 指数平滑法指在对过去数据加权平均时，对了离预测时间点越远的过去数据，赋予较小的权重，而越近的数据赋予越大的权重，过去数据对现在预测值的影响按照随距离呈指数形式衰减，其一般形式为 其中，为平滑系数。一般时间序列较稳定，取值可以小一些一般取；若时间序列数据起伏波动比较大则应取较大的值，一般。 拟合函数法通常采用最小二乘法对给定的形式的函数的参数进行估计，获得最终的解析函数表达式，从而可以计

3、算任何时间点的变量取值。3种最常用的趋势如下所示。 1、直线型趋势线；2指数型趋势线。3抛物线型趋势线 。 周期和季节的分析方法。 （1）自相关函数的分析方法对于一个时间序列，其自相关序列函数用于测定序列前后期数值之间的相关关系。 当序列中含有周期性成分时，其自相关序列也呈周期性变化，这一特性被用来判断时间序列中是否含有周期性变化成分以及周期的大小 matlab中提供了xcorr函数用来计算序列的自相关函数。也可以用autocorr函数来画系列的归一化后的自相关函数图。四、实验内容例一 分别利用3点平均法、3点滑动平均法、一次指数平滑法和直线拟合法对表1中的数据进行处理，重新计算在各时间点的估

4、计值 表一 某地区2001年2003年不同季度的住院人数年丰第一季度第二季度第三季度第四季度200125937234022020022754113532302003285428364243 clear;close all;clc;aa=259 372 340 220;275 411 353 230;285 428 364 243;bb=reshape(aa,1) %将数据按各年季度顺序排列为一列向量；figure,plot(bb,o-);title(某地区各季度住院人数（2001-2003年）);m1=bb %3点移动平均法for k=4:length(m1) %因为前3个点直接使用公式计算数

5、据不够，保持不变 m1(k)=mean(bb(k-3:k-1);endm2=bb %3点滑动平均法for k=2:length(m1)-1 %因为第一个点、最后一个点不能直接使用公式计算数据不够，保持不变 m2(k)=mean(bb(k-1:k+1);endm3=bb %指数平滑法alpha=0.3;for k=1:length(m1)-1 %因为第一个点、最后一个点不能直接使用公式计算数据不够，保持不变 m3(k+1)=alpha*bb(k)+(1-alpha)*m3(k);end%最小二乘直线拟合xx=1:length(bb)p=polyfit(xx,bb,1)m4=polyval(p,x

6、x)figure,plot(xx,m1,-.,xx,m2,-.*,xx,m3,-.d,xx,m4,-,xx,bb,-.);legend(三次移动平均,三次滑动平均,指数平滑法,直线拟合法,实测数据);bb = 259 372 340 220 275 411 353 230 285 428 364 243m1 = 259 372 340 220 275 411 353 230 285 428 364 243m2 = 259 372 340 220 275 411 353 230 285 428 364 243m3 = 259 372 340 220 275 411 353 230 285 428

7、 364 243xx = 1 2 3 4 5 6 7 8 9 10 11 12p = 1.7483 303.6364m4 = 305.3846 307.1329 308.8811 310.6294 312.3776 314.1259 315.8741 317.6224 319.3706 321.1189 322.8671 324.6154 图一 某地区各季度住院人数 图二 不同线性趋势预测方法对比 分析：（1）从上面我们就可以看到三次移动平均m1，三次滑动平均m2，指数滑动平均m3，直线拟和法m4等; （2）从图中可以看到，与原来的实测数据相比，几种方法都有效地减小了数据的波动性，更好地体现了

8、缓慢增加的变化趋势 不足之处是3中移动平均法明显存在着数据滞后，对现在的数据不敏感，指数平滑法较好地反应了趋势，直线拟合法从已测全部数据考虑 提取了最好的线性变化趋势。 例二 对表一中的数据进行处理，预测该地区2044年各季度的住院人数。clear;close all;clc; aa=259 372 340 220;275 411 353 230;285 428 364 243;bb=reshape(aa,1) %将数据按各年季度顺序排列为一列向量；figure,plot(bb,o-);title(某地区各季度住院人数（01-03年）);tt=1:length(bb)N=length(tt)f

9、igure,autocorr(bb);%序列的自相关图p=polyfit(tt,bb,1) %一次多项式拟合ytrend=polyval(p,tt) %线性项的预测值figure,plot(ytrend);yre=bb./ytrend; figure,plot(yre);yy=reshape(yre,4,)ysea=mean(yy); %求均值获得季节性指标yper=repmat(ysea,1,3); %12个季度的季节性指标ydesea=bb./yperp=polyfit(tt,ydesea,1) %一次多项式拟合tt2=1:16ytrend2=polyval(p,tt2); %线性项的预测

10、值yper=repmat(ysea,1,4); %16个季度的季节性指标ypredict=ytrend2.*yperfigure,plot(tt,bb,-.*,tt2,ypredict,-.o);legend(测量值,预测值)bb = 259 372 340 220 275 411 353 230 285 428 364 243tt = 1 2 3 4 5 6 7 8 9 10 11 12N = 12p = 1.7483 303.6364ytrend = 305.3846 307.1329 308.8811 310.6294 312.3776 314.1259 315.8741 317.622

11、4 319.3706 321.1189 322.8671 324.6154yy = 0.8481 0.8803 0.8924 1.2112 1.3084 1.3328 1.1007 1.1175 1.1274 0.7082 0.7241 0.7486ydesea = 296.4703 289.6869 304.8708 302.6207 314.7851 320.0573 316.5276 316.3762 326.2318 333.2957 326.3911 334.2583p = 3.7213 290.9428tt2 = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1

12、5 16ypredict = 257.4221 383.1700 336.9173 222.3316 270.4258 233.1527 283.4295 421.3991 370.1176 243.9739 296.4332 440.5136 386.7178 254.7950 图三 序列的自相关图 图4 由数据估计的线性趋势项 图5 去除线性趋势项后的数据 图6 某地区2001-2003年各季度住院人数的估计分析:（1）在本题中我们可以得到2004年四季度的住院人数的预测值分别为 296、440、386、254 (2）在这个题目中利用的额模型是乘法模型 （3）将不同年份相同的季度去除线性后

13、数据进行平均处理，即得到季节性指标。去除线性是 为了预测季节因素影响 （4）由序列的自相关图可以得到周期变化的步长为4,用一次多项式拟合可以估计为1.7485和b为303.6364即可以建立模型来描述该地区住院人数的变化趋势其中a、b是线性趋势项的参数，是季节性指标。 例3 脑电信号的分析 load eeg8s; %载入脑电数据testsig=eeg8s;figure,plot(testsig);xlim(1,1100)figure,histfit(testsig);muhat, sigmahat = normfit(testsig)figure,parcorr(testsig,50)%确定A

14、R模型阶数k=1;for p=5:20m=armax(testsig,p 0);FPEvalue(k)=m.EstimationInfo.FPE;pvalue(k)=p;k=k+1;endj=find(FPEvalue=min(FPEvalue);Porder=pvalue(j);figure,plot(pvalue,FPEvalue,-.o,pvalue(j),FPEvalue(j),*);xlabel(p-order)ylabel(FPE value)title(find the minimum FPE);m=armax(testsig,Porder,0);d1=m.a;p1=m.Nois

15、eVariance;H1,w1=freqz(sqrt(p1),d1);s = spectrum.periodogram;Hpsd = psd(s,testsig);%用周期图法估计信号的功率谱figure,plot(Hpsd); hold on;hp = plot(w1/pi,20*log10(2*abs(H1)/(2*pi),r); % Scale to make one-sided PSDset(hp,LineWidth,2);xlabel(Normalized frequency (times pi rad/sample)ylabel(One-sided PSD (dB/rad/samp

16、le)legend(PSD estimate of x,PSD of model output)m=armax(testsig,Porder,0);ye=predict(m,testsig,1);%模型预测输出figure, stem(testsig,ye1,1); xlabel(Sample time); ylabel(Signal value);legend(Original autoregressive signal,Signal estimate from linear predictor)p2 = norm(testsig-ye1,1),2)2/(length(testsig)-1)

17、;p1,p2muhat =0.0439sigmahat =1.3439ans =1.0158 0.9920 图7 单道脑电信号时间序列图 图8 脑电信号分布直方图图9 序列的PACF图图10 AR模型阶数图图11 脑电实测序列与AR模型的预测输出对比图 图12 脑电序列AR模型输出信号的功率谱分析：（1）载入单道脑电信号，首先用rand函数生成随机数，绘制其时序图；然后绘制直方图进行统计描述，看脑电信号是否服从正态分布，再进行正态分布拟合；判定AR模型的阶次。通过绘制序列的PACF图、再利用FPE准则确定AR模型阶数； 当模型阶数确定后，采用armax函数确定模型参数并进行预测。 （2）从图8

18、中可以看到脑电信号基本服从正态分布所以进行正态拟合的结果为均值=0.0439，方差=1.342 ，图9中可以看到的偏自相关函数约在20阶截尾。因利用FPE准则，对取值在5到20之间取值。图10中列出了不同阶次时的AR模型最终预测误差，从中可以看到P为11是FPE最小 ，因此确定AR模型的阶数为AR（11）。当模型确定后，可以采用armax函数确定模型参数并进行预测。模型预测的白噪声方差为0.9952 ，预测数据与实测数据间的方差为0.9853；在图12中可以看到模型预测的输出在频域上很好地拟合了实测的脑电数据。五、实验小结：（1）时间预测方法的基本思想是:预测一个现象的未来变化时，应该用现象的

19、过去行为来预测未来。（2）时间序列中的趋势项、周期项和季节变动项都可以看作确定性成分，除此外的部分则可以看作是随机的成分，时间分析方法可以分为两类：确定型时间序列分析方法和随机时间序列分析方法。（3）时间序列分析首先要确定序列的统计特性是否平稳从而选择不同的分析方法(4)时间序列分析的目的是1、描述事物在过去发生的规律，2分析事物发展的规律 3、对事物的发展规律进行预测或加于控制（5）时间序列通常可以用加法和乘法两种模型来表示。 (6)对于确定型的时间序列，如何看待其中不确定的成分？ 不确定成分是噪音，与确定成分没有任何相关性（7）对于随机时间序列，如何看待其中的确定成分? 随机时间序列中的确定成分是含有一定相关结构的白噪声线性组合。（8）对于平稳序列，简述进行AR建模的步骤. 首先确定AR的模型的阶次，然后确定AR模型的参数1. (9)时间序列分析是根据其时间历程是否具有规律性重复出现的周期成分分为周期性和非周期性的。 手写签名：