时间序列分析 2.docx
《时间序列分析 2.docx》由会员分享，可在线阅读，更多相关《时间序列分析 2.docx（21页珍藏版）》请在冰豆网上搜索。
时间序列分析2
评分
2014—2015学年度第3学期
一、实验目的
在生物医学领域，时间序列分析被广泛用于生理信号分析（如心电信号、脑电信号）流行病学研究（如传染病发病情况的预测）医学图像特征提取、（如fMEI图像序列中功能活动的激活点提取）等方面。
本次实验要求熟悉时间系列分析的移动平均法、滑动平均法、自相关函数分析法，了解脑电信号分析的matlab示例。
二、实验环境
1、硬件配置：
处理器：
Intel（R）Pentium（R）CPU@1.80GHz3.1880GHz
安装内存（RAM）：
4.00GB系统类型：
64位操作系统
2、软件环境：
MATLABR2010b软件
三、实验内容
（包括本实验要完成的实验问题及需要的相关知识简单概述）
1.实验原理
趋势分析方法
（1）移动平均法。
设某一时间序列为
，则
时刻的预测值为
（2）滑动平均法
式中；为t点的滑动平均值；
为单侧平滑时距
（3）指数平滑法
指数平滑法指在对过去数据加权平均时，对了离预测时间点越远的过去数据，赋予较小的权重，而越近的数据赋予越大的权重，过去数据对现在预测值的影响按照随距离呈指数形式衰减，其一般形式为
其中，
为平滑系数。
一般时间序列较稳定，
取值可以小一些一般取
；若时间序列数据起伏波动比较大则
应取较大的值，一般
。
拟合函数法通常采用最小二乘法对给定的形式的函数的参数进行估计，获得最终的解析函数表达式，从而可以计算任何时间点的变量取值。
3种最常用的趋势如下所示。
1、直线型趋势线
；2指数型趋势线
。
3抛物线型趋势线
。
周期和季节的分析方法。
（1）自相关函数的分析方法
对于一个时间序列，其自相关序列函数用于测定序列前后期数值之间的相关关系。
当序列中含有周期性成分时，其自相关序列也呈周期性变化，这一特性被用来判断时间序列中是否含有周期性变化成分以及周期的大小
matlab中提供了xcorr函数用来计算序列的自相关函数。
也可以用autocorr函数来画系列的归一化后的自相关函数图。
四、实验内容
例一分别利用3点平均法、3点滑动平均法、一次指数平滑法和直线拟合法对表1中的数据进行处理，重新计算在各时间点的估计值
表一某地区2001年—2003年不同季度的住院人数
年丰
第一季度
第二季度
第三季度
第四季度
2001
259
372
340
220
2002
275
411
353
230
2003
285
428
364
243
clear;
closeall;
clc;
aa=[259372340220;275411353230;285428364243];
bb=reshape（aa',[],1）%将数据按各年季度顺序排列为一列向量；
figure,plot（bb,'o-'）;
title（'某地区各季度住院人数（2001-2003年）'）;
m1=bb%3点移动平均法
fork=4:
length（m1）%因为前3个点直接使用公式计算数据不够，保持不变
m1（k）=mean（bb（k-3:
k-1））;
end
m2=bb%3点滑动平均法
fork=2:
length（m1）-1%因为第一个点、最后一个点不能直接使用公式计算数据不够，保持不变
m2（k）=mean（bb（k-1:
k+1））;
end
m3=bb%指数平滑法
alpha=0.3;
fork=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,xx）
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
243
m1=
259
372
340
220
275
411
353
230
285
428
364
243
m2=
259
372
340
220
275
411
353
230
285
428
364
243
m3=
259
372
340
220
275
411
353
230
285
428
364
243
xx=
1
2
3
4
5
6
7
8
9
10
11
12
p=
1.7483303.6364
m4=
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）从图中可以看到，与原来的实测数据相比，几种方法都有效地减小了数据的波动性，更好地体现了缓慢增加的变化趋势不足之处是3中移动平均法明显存在着数据滞后，对现在的数据不敏感，指数平滑法较好地反应了趋势，直线拟合法从已测全部数据考虑提取了最好的线性变化趋势。
例二对表一中的数据进行处理，预测该地区2044年各季度的住院人数。
clear;
closeall;
clc;
aa=[259372340220;275411353230;285428364243];
bb=reshape（aa',[],1）%将数据按各年季度顺序排列为一列向量；
figure,plot（bb,'o-'）;
title（'某地区各季度住院人数（01-03年）'）;
tt=[1:
length（bb）]'
N=length（tt）
figure,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./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（'测量值','预测值'）
bb=
259
372
340
220
275
411
353
230
285
428
364
243
tt=
1
2
3
4
5
6
7
8
9
10
11
12
N=
12
p=
1.7483303.6364
ytrend=
305.3846
307.1329
308.8811
310.6294
312.3776
314.1259
315.8741
317.6224
319.3706
321.1189
322.8671
324.6154
yy=
0.84810.88030.8924
1.21121.30841.3328
1.10071.11751.1274
0.70820.72410.7486
ydesea=
296.4703
289.6869
304.8708
302.6207
314.7851
320.0573
316.5276
316.3762
326.2318
333.2957
326.3911
334.2583
p=
3.7213290.9428
tt2=
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
ypredict=
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）将不同年份相同的季度去除线性后数据进行平均处理，即得到季节性指标。
去除线性是
为了预测季节因素影响
（4）由序列的自相关图可以得到周期变化的步长为4,用一次多项式拟合可以估计
为1.7485
和b为303.6364即可以建立模型
来描述该地区住院人数的变化趋势其中a、b是线性趋势项的参数，
是季节性指标。
例3脑电信号的分析
>>loadeeg8s;%%%载入脑电数据
testsig=eeg8s;
figure,plot（testsig）;xlim（[1,1100]）
figure,histfit（testsig）;
[muhat,sigmahat]=normfit（testsig）
figure,parcorr（testsig,50）
%%%确定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）,'*'）;
xlabel（'p-order'）
ylabel（'FPEvalue'）
title（'findtheminimumFPE'）;
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'）
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）;
[p1,p2]
muhat=
0.0439
sigmahat=
1.3439
ans=
1.01580.9920
图7单道脑电信号时间序列图
图8脑电信号分布直方图
图9序列的PACF图
图10AR模型阶数图
图11脑电实测序列与AR模型的预测输出对比图
图12脑电序列AR模型输出信号的功率谱
分析：
（1）载入单道脑电信号，首先用rand函数生成随机数，绘制其时序图；
然后绘制直方图进行统计描述，看脑电信号是否服从正态分布，再进行正态分布拟合；判定AR模型的阶次。
通过绘制序列的PACF图、再利用FPE准则确定AR模型阶数；当模型阶数确定后，采用armax函数确定模型参数并进行预测。
（2）从图8中可以看到脑电信号基本服从正态分布所以进行正态拟合的结果为均值=0.0439，方差=1.34^2，图9中可以看到的偏自相关函数约在20阶截尾。
因利用FPE准则，对取值在5到20之间取值。
图10中列出了不同阶次时的AR模型最终预测误差，从中可以看到P为11是FPE最小，因此确定AR模型的阶数为AR（11）。
当模型确定后，可以采用armax函数确定模型参数并进行预测。
模型预测的白噪声方差为0.9952，预测数据与实测数据间的方差为0.9853；在图12中可以看到模型预测的输出在频域上很好地拟合了实测的脑电数据。
五、实验小结：
（1）时间预测方法的基本思想是:
预测一个现象的未来变化时，应该用现象的过去行为来预测未来。
（2）时间序列中的趋势项、周期项和季节变动项都可以看作确定性成分，除此外的部分则可以看作是随机的成分，时间分析方法可以分为两类：
确定型时间序列分析方法和随机时间序列分析方法。
（3）时间序列分析首先要确定序列的统计特性是否平稳从而选择不同的分析方法
（4）时间序列分析的目的是1、描述事物在过去发生的规律，2分析事物发展的规律3、对事物的发展规律进行预测或加于控制
（5）时间序列通常可以用加法和乘法两种模型来表示。
（6）对于确定型的时间序列，如何看待其中不确定的成分？
不确定成分是噪音，与确定成分没有任何相关性
（7）对于随机时间序列，如何看待其中的确定成分?
随机时间序列中的确定成分是含有一定相关结构的白噪声线性组合。
（8）对于平稳序列，简述进行AR建模的步骤.
首先确定AR的模型的阶次，然后确定AR模型的参数
1.（9）时间序列分析是根据其时间历程是否具有规律性重复出现的周期成分分为周期性和非周期性的。
手写签名：
升级会员 