时间序列分析.docx
《时间序列分析.docx》由会员分享，可在线阅读，更多相关《时间序列分析.docx（13页珍藏版）》请在冰点文库上搜索。
时间序列分析
一、
实验目的
在生物医学领域，时间序列分析被广泛用于生理信号分析（如心电信号、脑电信号的研究）、流行病学研究（如传染病发病情况的预测）、医学图像特征提取（如fMRI图像序列中功能活动是脑内的激活点提取）等方面。
本次试验要求熟悉掌握时间序列分析的移动平均法、滑动平均法、自相关函数分析法，了解脑电信号分析的Matlab示例。
二、实验环境
1、硬件配置：
处理器（Intel（R）Core（TM）i7-3770CPU@3.40GHz3.40GHz、CD-ROM驱动器、鼠标、内存1GB（1024MB）、64位操作系统
2、软件环境：
MATLABR2010b
三、实验内容
（包括本实验要完成的实验问题及需要的相关知识简单概述）
1.试验原理
平滑方法属于比较古典的方法，主要有移动平均法、滑动平均法和指数平滑法，基本思想是取待预测点前（后）几个值得平均或加权平均作为预测值，往往在不易判别趋势项函数形式的情况下使用。
（1）移动平均法。
设某一时间序列为y1，y2，…,yi,则t+1时刻的预测值为
式中：
为t点的移动平均值；n称为移动时距。
（2）滑动平均法
式中：
为t点的滑动平均值；l为单侧平滑时距。
（3）自相关函数分析法
对于一个时间序列，其自相关函数用于测定序列前后期数值之间的相关关系。
当序列中含有周期性成分时，其自相关序列也呈周期性变化，这一特性常被用于判断时间序列中含有周期性变化的成分以及周期的大小。
MATLAB中提供了xcorr函数来计算序列的自相关函数。
也可以用autocorr函数来画序列的归一后的自相关函数图。
2.实验内容
（1）P132，例5.1
分别利用3点移动平均法、3点滑动平均法、一次指数平滑法和直线拟合法对表1中的数据进行处理，重新计算在各时间点的估计值。
（2）P135，例5.2
对表1中的数据进行处理，预测该地区2004年各季度的住院人数。
（3）P145，例5.5
脑电信号的分析
表1某地区2001年-2003年不同季节的住院人数
年份
第一季度
第二季度
第三季度
第四季度
2001
259
372
340
220
2002
275
411
353
230
2003
285
428
364
243
四、实验结果与分析
%%%exam51.m
clear;
closeall;
clc;
aa=[259372340220;275411353230;285428364243];
bb=reshape（aa',[],1）;%%将数据按各年季度顺序排列为一列向量；
figure,plot（bb,'o-'）;title（'某地区各季度住院人数（01－03年）'）;%%绘制某地区不同季度的住院人数如图1所示
%%%三点移动平均法
m1=bb;
fork=4:
length（m1）%%因为前三个点直接使用公式计算数据不够，保持不变
m1（k）=mean（bb（k-3:
k-1））;
end
%%%三点滑动平均法
m2=bb;
fork=2:
length（m1）-1%%因为第1个点、最后一个点不能直接使用公式计算数据不够，保持不变
m2（k）=mean（bb（k-1:
k+1））;
end
%%%指数平滑法
m3=bb;
alpha=0.3;
fork=1:
length（m1）-1%%因为第1个点、最后一个点不能直接使用公式计算数据不够，保持不变
m3（k+1）=alpha*bb（k）+（1-alpha）*m3（k）;
end
%%%最小二乘直线拟和
xx=[1:
length（bb）]';
p=polyfit（xx,bb,1）;
m4=polyval（p,xx）;
figure,plot（xx,m1,'-.',xx,m2,'-.*',xx,m3,'-.d',xx,m4,'-',xx,bb,'-.^'）;
legend（'三次移动平均','三次滑动平均','指数平滑法','直线拟和法','实测数据'）;%%绘制不同线性趋势预测方法对比如图2所示
图1图2
结果分析：
首先这一题的实验数据的选取要根据实验原理中的各式子中自变量的取值范围来选取；另：
据图2示，与原来的实测数据相比，四种分析方法都有效地减小了数据的波动性，更好地体现了缓慢增长的变化趋势，但3次移动平均明显存在着数据滞后，对现在时刻的数据不敏感，而指数平滑法较好的反映了趋势，而直线拟合法从已测的全部数据考虑，提取最佳的线性变化趋势。
三点移动平均法在各时间点估计值：
259372340323；310278302346；331289314359
三点滑动平均法在各时间点估计值：
259323310278；302346331289；314359345243
指数平滑法在各时间点估计值：
259259292307；280279318328；299295334343
最小二乘直线拟和在各时间点估计值：
305307308310；312314315317；319321322324
%%%exam52.m
clear;
closeall;
clc;
aa=[259372340220;275411353230;285428364243];
bb=reshape（aa',[],1）;%%将数据按各年季度顺序排列为一列向量；
figure,plot（bb,'o-'）;title（'某地区各季度住院人数（01－03年）'）;%%绘制图1某地区不同季度的住院人数
tt=[1:
length（bb）]';
N=length（tt）;
figure,autocorr（bb）;%%绘制序列的自相关图如图3所示
图3
p=polyfit（tt,bb,1）;%%一次多项式拟合
ytrend=polyval（p,tt）;%%线性项的预测值
figure,plot（ytrend）;
yre=bb./ytrend;
figure,plot（yre）;%%绘制数据中的线性趋势项提取过程图如图4所示
yy=reshape（yre,4,[]）;
ysea=mean（yy'）;%%求均值获得季节性指标
yper=repmat（ysea,1,3）;%%12个季度的季节性指标
ydesea=bb./yper';
p=polyfit（tt,ydesea,1）;%%一次多项式拟合
tt2=[1:
16]';
ytrend2=polyval（p,tt2）;%%线性项的预测值
yper=repmat（ysea,1,4）;%%16个季度的季节性指标
ypredict=ytrend2.*yper';
figure,plot（tt,bb,'-.*',tt2,ypredict,'-.o'）;legend（'测量值','预测值'）%%绘制某地区2001-2004年各季度住院人数估计图如图5所示
（a）由线性项预测的数据变化趋势（b）消除线性趋势项后的数据图
图4
图5
分析：
由图1可知该地区各季度住院人数随年份缓慢上升，并且具有季节因素；如图3所示的序列的自相关图可以验证各季度住院人数有周期变化，步长为4个采样步长，之后建立函数模型描述该地区住院人数的变化趋势，首先提取线性趋势项来估计线性趋势项参数a和b，将不同年份相同季度的去除线性项后数据进行平均处理，得到季节性指标，用原始数据除以季节性指标后重新估计线性趋势项。
最后进行预测时，用线性项预测值乘以季节性指标得到各时间点预测值。
2001-2004年各季度住院人数预测值：
257383336222；%%第一季度
270402353233；%%第二季度
283421370243；%%第三季度
296440386254；%%第四季度
%%%exam55.m
clear;
closeall;
clc;
loadeeg8s;%%%载入脑电数据
testsig=eeg8s;
figure,plot（testsig）;xlim（[1,1100]）%%绘制单道脑信号时间序列图如图6所示
figure,histfit（testsig）;%%绘制脑信号分布直方图如图7所示
图6图7
%%图6是依据单道脑电信号绘制的时序图；图7是脑电信号的直方图，且从直方图可看出脑电信号基本服从正态分布。
[muhat,sigmahat]=normfit（testsig）
figure,parcorr（testsig,50）%%绘制序列的PACF图如图8所示
图8
%%%确定AR模型阶数
k=1;
forp=5:
20
m=armax（testsig,[p0]）;
FPEvalue（k）=m.EstimationInfo.FPE;
pvalue（k）=p;
k=k+1;
end
j=find（FPEvalue==min（FPEvalue））;
Porder=pvalue（j）;
figure,plot（pvalue,FPEvalue,'-.o',pvalue（j）,FPEvalue（j）,'*'）;%%绘制利用FPE准则确定AR模型阶数图如图9所示
xlabel（'p-order'）
ylabel（'FPEvalue'）
title（'findtheminimumFPE'）;
图9
m=armax（testsig,[Porder,0]）;
d1=m.a;
p1=m.NoiseVariance;
[H1,w1]=freqz（sqrt（p1）,d1）;
s=spectrum.periodogram;
Hpsd=psd（s,testsig）;%%用周期图法估计信号的功率谱
figure,
plot（Hpsd）;
holdon;
hp=plot（w1/pi,20*log10（2*abs（H1）/（2*pi））,'r'）;%Scaletomakeone-sidedPSD
set（hp,'LineWidth',2）;
xlabel（'Normalizedfrequency（\times\pirad/sample）'）
ylabel（'One-sidedPSD（dB/rad/sample）'）
legend（'PSDestimateofx','PSDofmodeloutput'）%%绘制脑电序列及AR模型输出信号功率谱图如图10所示
图10
m=armax（testsig,[Porder,0]）;
ye=predict（m,testsig,1）;%%模型预测输出
figure,
stem（[testsig,ye{1,1}]）;
xlabel（'Sampletime'）;
ylabel（'Signalvalue'）;
legend（'Originalautoregressivesignal','Signalestimatefromlinearpredictor'）
p2=norm（（testsig-ye{1,1}）,2）^2/（length（testsig）-1）;%%绘制脑电实测序列与AR模型的预测输出图如图11所示
图11
[p1,p2]
%%实验结果输出
muhat=
0.0439
sigmahat=
1.3439
ans=
1.01580.9920
分析：
模型预测的白噪声方差为1.0158，而预测数据与实测数据间残差的方差为0.9920。
模型预测值与实际测量值如图10所示。
由于信号的统计特性，从时域图上很难看出模型的预测效果，因此绘制脑电信号的功率谱与AR系统预测信号的功率谱进行对比，如图11所示，可以看到模型预测的输出在频域上很好地拟合了实测的脑电数据。
载入单道脑电信号，首先绘制其时序图，如图6所示。
然后绘制直方图进行统计描述，如图7所示。
从直方图来看脑电信号基本服从正态分布，因此可进行正态分布的参数估计。
进行正态分布拟合的结果为：
均值=0.0439，方差=1.34^2。
之后用高阶的AR模型进行逼近，首先应判定AR模型的阶次。
图8为序列的PACF图，从图中可以看到，序列的偏自相关函数约在20阶后截尾。
因此利用FPE准则，对P取值在5-20之间进行判断。
从图9中看到P为11时FPE最小，因此确定AR模型的阶数为AR（11）。
模型阶数确定后，采用armax函数确定模型参数并进行预测。
五、实验小结：
本次试验做的是时间序列分析：
1.时间序列分析是以预测一个现象的未来变化时，应该用该现象的过去行为来预测未来。
也就是通过时间序列的历史数据探索现象变化的规律，并将这种规律衍生出对未来的预测。
2.时间序列分析一般分为：
（1）序列的预处理
（2）描述性时序分析（3）统计时序分析（4）判断时间序列的组合成分，对不同成分采取不同的分析方法（5）利用不同成分的模型进行预测后叠加，得到最终的时间序列分析预测值。
3.时间序列分析是根据其时间历程是否具有规律性重复出现的周期成分分为周期性和非周期性的。
4.对于随机型的时间序列，首先要判断序列的统计特性是否平稳性。
A．对于平稳性时间序列的分析，可以在时域进行，也可以在频域进行。
即功率谱估计（周期图法、平均周期图法、平滑周期图法和平均平滑图法）和现代建模法（AR模型、MA模型、ARMA模型）。
B．对于非平稳性时间序列，本次实验主要涉及ARIMA模型，即采用差分的方法来做平稳化处理，使非平稳序列变成平稳序列，然后按照平稳序列建立的ARMA模型进行预测，最后再通过做差分的逆操作得到原序列的预测值。
5.动平均、滑动平均、指数平滑中对于数据的选择是一个难点，但都可以根据各自的运算公式发现规律。
6.其中例题5.5是一个较为完整的对非平稳性时间序列的分析示例包括对序列进行画图描绘、统计分析、确定AR模型的阶数和参数、确定AR模型的预测输出，由于在时域不易判断所采用的模型是否有效地描述了随机过程，因此采用功率谱来对照分析模型的预测结果。
手写签名：
升级会员 