文稿系统辨识

1、大连理工大学系统辨识课题目名称：Generalized additive sforlocationand(GAMLSSin评学大连理工大学系统辨识课题目名称：Generalized additive sforlocationand(GAMLSSin评学院（系电专业仪器仪表学位类别学术型 专业型 姓名号： 导师秦完成日期: 2015年6月81 研究背景及意时间序列(timeseries)1.1 1 研究背景及意时间序列(timeseries)1.1 图地磁扰动强度Dst指数时间序时间序列分析(Time1.2 描建决策和控1.2 1.2 描述：用图形(1.1)来描述时间序列，或应用样本的自协方差、

2、样本相关即是指从Dst阳风速度和等离子密度等变量呈现了整数型泊松分布或实数型系统中， 因此，开展非平稳时间序列的研究具有重要的理论意义和应用价值产生时间序列的系统看作随机系统，并将所获得的数据集 =| = 0,1,2,的观测值量1.2 描述：用图形(1.1)来描述时间序列，或应用样本的自协方差、样本相关即是指从Dst阳风速度和等离子密度等变量呈现了整数型泊松分布或实数型系统中， 因此，开展非平稳时间序列的研究具有重要的理论意义和应用价值产生时间序列的系统看作随机系统，并将所获得的数据集 =| = 0,1,2,的观测值量就是要根据观测到的数据集建立关于的条件概率密度模型(|AI,) 。其中，AI

3、于数据集国内外研究及发展现1.1表JenkinsetEngleand Timeysis：Forecastingand(AutoRegressiveConditional TimeSeriesandysisPanditand安、 NelderandMcCullagh(1989)表JenkinsetEngleand Timeysis：Forecastingand(AutoRegressiveConditional TimeSeriesandysisPanditand安、 NelderandMcCullagh(1989)6杨叔子AlcazarandSharman(1996)9etGeneralize

4、dLinear。Rosipalnet。GAMLSS(GeneralizedadditivesRigbyetlocationscaleand)Grangerand图通过图 1.3 AR ARMA(AutoRegressive and Moving Average ARARX(AutoRegressive exogenous) ARMA 通过图 1.3 AR ARMA(AutoRegressive and Moving Average ARARX(AutoRegressive exogenous) ARMA ARXMA(AutoRegressiveexogenousandMovingAverage

5、况下，图1.3中的1.4 当状态变量不可观测时通常用图1.4所示的状态空间方程来描述时间序列中，为独立同分布的。虽然如HammersteinWiener模型18AR 类模型中要求数据的方差为稳态的这一问题，Robert Engle 归条件异方差(AutoRegressiveConditionalHeteroskedasticity2 = 0 + 其中有效的解决一些实际问题。ARCH 2003 经济学奖的计量经济学成果之一。但由于其规定了此在理论与实际应用上仍具有一定的局限性262 = 0 + 其中有效的解决一些实际问题。ARCH 2003 经济学奖的计量经济学成果之一。但由于其规定了此在理论与

6、实际应用上仍具有一定的局限性26近任意的分布15。在信号处理领域，Yoonn 出了类似的方案16。由于核密度估计模型充分考虑了随机输出变量概率分布的非对称 较大。为解决这一问题，本研究将基于广义线性模型(Generalized linear ，GLM模型17 。在统计学上，GLM 模型是一种受到广泛应用的线性回归模式。GLM 模型是基于指 logistic 回归以及泊松回归等经典的回归模型。在某些条件下，GLM 究非平稳的时间序列数据进行建模。AIC(Akaike information criterion)19Information Criterions)201.2所示。此后，Konishi

7、GIC24，GIC 外，Lasso Zou 等还将1-及2-Elastic net 当的正则化参数，Lasso 度以及信息量准则问题也受到一定的关注。除了1-外，Lasso Zou 等还将1-及2-Elastic net 当的正则化参数，Lasso 度以及信息量准则问题也受到一定的关注。除了1-等1/2-Bridge 表AIC 、BIC对(Goodness of Fit)的一本文主要内容及创图1.5 计出的模型利用决定系数进行评价，并利用BC信息量准则进行估计，选用BIC的原因将在第三章给出；最后，将研究方法用于空间天气预报，并对所提方法进行反馈修正。图1.5 计出的模型利用决定系数进行评价，

8、并利用BC信息量准则进行估计，选用BIC的原因将在第三章给出；最后，将研究方法用于空间天气预报，并对所提方法进行反馈修正。本研究的主要创新在于将基于GLM所模型可以对决定概率分布性质的数学期望，方差，偏度(Skewness)与峰度分布，而根据本研究的结果发现，使用基于2 回归意义下的GAMLSS 模GLM 时的情况，对于像服从于分布的时间序列，有数学期望，方差度三GLM GLM 、Akantziliotou 、Stasinopoulos 2002 年提 additive s)28的2 回归意义下的GAMLSS 模GLM 时的情况，对于像服从于分布的时间序列，有数学期望，方差度三GLM GLM

9、、Akantziliotou 、Stasinopoulos 2002 年提 additive s)28的局限性后GAMLSS模型27(2.1)其中 = (1233) = () () = = +1 11 1 () = = + 2 22 2 () = = + 3 33 3 () = = + 4 44 4 经过研究 发现，GAMLSS是一种(半独立式)参变的回归模型，这意味着假设响GAMLSS 中对于响应变量要服从指数族分布这一要求是很变量或是随机影响都是可以的。因此，GAMLSS 特别适合那些 GAMLSS模型中，如果令 = , = ( , , , ) ， 1211)，1是= ( )的逆如果 是单

10、数的话， 将服从与exp( 成 2比例的错误的先验密度函数。由公式(2.1)GAMLSS GAMLSS () ( ) = = + ( 其中 = 1,2,3,4时分别为, , , ，()为关于解释变量未知函数，为在GAMLSS()1() = = 如下式所示，模型(2.5)可以扩展为具有 ( ) = = ( , ) + 如下式所示，模型(2.5)可以扩展为具有 ( ) = = ( , ) + ( (3.26)称为非线性半参数累加模型。如果 = 可以得到非线性参数GAMLSS() = = (,如果() = (2.7)与(2.5)GAMLSS (3.27)中的() GAMLSS 模型可以看成线性和非线

11、性参数模型的结合。因此 称(2.7)或(2.5)GAMLSS 模型。在GAMLSS 模型中，本研究利用罚似然函数(2.8)来进行模型参数与 = 1 2其中=log (|)。在本研究GAMLSS3 进一步研究设3.1 经典时间序列模表3.13.1 AR可以噪声较为复杂的ARMA在时间序列分析中很多模型都可以的用状态空交替计算 = 1(1 )和 = (3 进一步研究设3.1 经典时间序列模表3.13.1 AR可以噪声较为复杂的ARMA在时间序列分析中很多模型都可以的用状态空交替计算 = 1(1 )和 = (滤波多时间序列分析都可以转换化为状态空间模3.2 GLM模独立。在数据处理过程中认为响应的过

12、去值的条件期望 = E|1是一个关于协变为响应的过去值的条件期望 于跃迁概率为 的二元时间序列来说， 对协变量进行线性回归分析 = 的Poisson或者Gamma 模型6(3.1)，这是一种严格服从系统性和随机1() = 2() = 这里 = 1, ( ( ; ,) = exp+ ( ;这里 = 1, ( ( ; ,) = exp+ ( ;其中 () ， 离散参数是已知的权值或者过去的权值， 是分布的自然2 系统性部分：对于 = 1,有单调函数g(.() = = = 函数g(. )是当变量向量 接函数。当一维响应向量 = 0 + 11 + 22 + 3cos 或 = 0 + 11 +22 +

13、31 + = + )+ ) 其中, 1 = ,1(1),(),1(1), = (,1,1,() = = + ) ) + 其中= () ，因为包含了自回归和滑动平均两个部分，被称为GLARMA(Generalized Linear Autoregressive Moving Average)GARMA(Generalized Autoregressive Moving Average)。() = = 目前Box-Cox = 1 =0 表示对数线性模型。需要注意的是当未知参数为相乘关系时，在估计GARMA模型时在上述分析中，关于1， = E|1 的解释如下。假设= ( )序列() = = 目前Bo

14、x-Cox = 1 =0 表示对数线性模型。需要注意的是当未知参数为相乘关系时，在估计GARMA模型时在上述分析中，关于1， = E|1 的解释如下。假设= ( )序列p 维向量 ( = 1,)。由1，1确定的11 = (1,2 ,1,2 知道对于 = 1,有(|1 = 1同时对： ( (; ,) = = 将等式(3.8)两边同时对：(;,|1) (;,|1) (; , |1) = = ( ( (;,|1) = (;,|1)(; , |1) = () (; , |1) =：Var|1 = 因为Var|1 0，是单调的，因此可以将(3.10) = 因此由连接函数(因为Var|1 0，是单调的，因

15、此可以将(3.10) = 因此由连接函数( ，1 ()1可以导出 = = 通过以上分析它的特点如表3.1表在统计学上，GLM此模型假设实验者所量GLM测的随量的分布函GLM数与实验中系统性效应 (即非随机的效应)可由一 连 接 函 数 (link function)建立起可以解GLMGLM 模型服从条件概率密度函数为(|21AI)的概率分布。其中AI为辅助信息量，然而在大多数情况下 AI AI ) = = (|1,2,1)(,1,2,= (|1,2,1)(1|1,2,2)(,1,2,在有协变量参与时似然函数() = (1122) 其中 = (112 其中 = (112211)， = (1122

16、11)() = (; , ()=log(;, =+( ; ( ) ( = +(;其中() ()1 = 11()，并且从(3.9)可以推出 = ( ,1 = 其中 = 1,。由(3.9)和(3.10) = 1( Var 又因为= = = Var|1其中 = 1,。基于以上分析定义() = = Var|1其中 = 1,。基于以上分析定义() () () = =1 其中2() = Var|1。对于 = 1,定义() () = =1 2() E = 1 2() 假设E() = 0，当 () () = E() = log() = 由此可以解出AR 模型中，通常定义其噪声ARGAMLSSDst数据进行模型

17、估计与选择时，首先假设其服分布利用AR 模型进行数据的建模3.1 为方差2分布概率密度函数曲线。如果(量1() ( 由此可以看分布是典型的指数族分布。由图3.1 可以看曲线关于轴对称，并且|3.1 为方差2分布概率密度函数曲线。如果(量1() ( 由此可以看分布是典型的指数族分布。由图3.1 可以看曲线关于轴对称，并且|的值越大，()的值越趋近于零，当| 时() 0 AR 模型中定义噪声(0,2)望和方差，因此在进行模型的估计与选择是只需估计和。GAMLSS = +1 log()=+0基于AR由于GAMLSS服从分布时进行对数学期望，标准差度称设(0,1)，2且和 度为 的 分布29。 其中，

18、 如果是 变量， 其分布称为12 i.i.d.(0,1) =度为的2变量，其分布称度为的2分布，记为2 度为 的 分布29。 其中， 如果是 变量， 其分布称为12 i.i.d.(0,1) =度为的2变量，其分布称度为的2分布，记为2量，其概率密度函数可以写为(3.21)3.2 2) 2(1n3.2 3.1 分布曲线很相似。二者均为关于原点对称，单峰值的偶函数，并在 = 0但分布相比需要注意的是，()当且仅当 1) + )22n1() ( )( 当 2时， + )22n1() ( )( 当 2时，() = 0；当 3时，() 当 = 1时，t11() = (1 + 其中 ()当 时，变量的极限

19、分布为(0,1)Dst GAMLSS = +1 log( ) = +2 0 log() = +3 0模型的评价与选determination)(3.22) 2 = 1 决定系数2是相关系数的平方， 数据对随机模型的拟合程度的一个随机变量。决定系数通常基于相关信息对输出值进行 或者对 进行检验。目前根据不在以上两种情况中，2 (0,1)。当性回归模型，为解决这一问题，Cox 和 Snell 在，()上节给出了当自回归阶数给定情况下模型参数数AIC BIC 来确定最优的自回归阶数(延时) = 1, BICBIC 最小时的4由上一章的内不难发现，本研究是将GLM 模型扩展到时间序列对非 风与Dst指

20、数的相上节给出了当自回归阶数给定情况下模型参数数AIC BIC 来确定最优的自回归阶数(延时) = 1, BICBIC 最小时的4由上一章的内不难发现，本研究是将GLM 模型扩展到时间序列对非 风与Dst指数的相风与Dst指数介Dst(Disturbance Storm Time)Worldenter eomagnetism, Kyoto 2012的4.1 的左侧，右侧为4.1 Nagelkerke 在1991 年对其有了进一步的研究25，由于关于非线性模型22 = 1 ()11 风的速度一般为200800km/s极风一词是由尤金派克301950 1960 于等、4.2Dst4.2 Dst D

21、st工Dst仿真实验结11 风的速度一般为200800km/s极风一词是由尤金派克301950 1960 于等、4.2Dst4.2 Dst Dst工Dst仿真实验结4.3DstDstDst指数数据的曲4.3 4.4 为数据的直方图。4.4Dst4.4Dst由于 分布的似然函数较为简单，并且容易计算。因此本研究首先考虑数据服从4.4 100 的分析了解到分布有期望、方差、 度这三个对数据分布情况影响较大的统计量。在本研究中尝试对期望、方差、 度分别进行动态估计，其具体过程如下：假设数据的方差度为定常，对期望进行动态建模假设数据度为定常，分别对期望、方差进行动态建模假设数据的期望、方差度均为时变的

22、，分别对三者进行动态建模本研究中利用R中的GAMLSS数据包计算地磁Dst指数时间序列的模型参数、决定系数以及BIC 选择结果。方差为定常分布数据的仿真实本仿真实验中，令最大延时 = 504.5AR4.5Estimate 4.5AR4.5Estimate 为模型参数的估计值，Std.Error SBCBICBIC值。通过仿真结果可以看出，BIC选择出的最优的最大延时 = 34.6 为模型拟合后的曲线。图AR模型拟合曲4.7 图时变AR模型的仿真结BIC 信息量准则选择的最优的最大延时 = 图时变AR模型的仿真结BIC 信息量准则选择的最优的最大延时 = 4.2.1 BIC 4.8 图时变的AR

23、模型拟合曲方差为常数的服从于4.9 4.9 分布仿真结果(方差定常于4.9 4.9 分布仿真结果(方差定常 = 504.9度都是常数的情况下，BIC4.10 4.10 分布拟合曲线(方差定常度为定常的服从。由4.2.3 节。由4.2.3 节度都是常数，通过BIC 选择出的图4.11 为估计结果。图4.11 分布的仿真结果度定常BIC对数据的方差进行估计时，BIC 信息量准则的选择结果 = 2，可以认为这是一个较为4.12 为拟合曲线。4.12 服从在及 = 504.13 时变的在及 = 504.13 时变的4.13 BIC选择出的最优延时为3BIC择结果为 = 2BIC 选择出的最优延时为 =

24、 24.14 仿真小BIC4.1 3.4.1 BICAR33AR32仿真小BIC4.1 3.4.1 BICAR33AR323224.1 AR 对于服从分布的模型，只有期望时变的这种情况，BIC BICBox G E P, Jenkins G M, Reinsel G C. Time &Sons,ysis: forecasting and controlM. John Engle R F. Autoregressive conditional heteroscedasticity with estimates of the variance of KingdominflationJ.Econom

25、etrica:JournaloftheEconometricSociety,1982:987- Box G E P, Jenkins G M, Reinsel G C. Time &Sons,ysis: forecasting and controlM. John Engle R F. Autoregressive conditional heteroscedasticity with estimates of the variance of KingdominflationJ.Econometrica:JournaloftheEconometricSociety,1982:987- Pand

26、itSM,WuSM.TimeseriesandysispplicationsM.NewYork:Wiley,安，杜金观，潘一民时间序列的分析与应用Bollerslev T. Generalized autoregressive conditional heteroskedasticityJ. Journal of 1986,31(3):307-McCullagh P. Generalized linear sJ. European Journal of Operational Research, 1984, BaughmanDR.Neuralnetworksinsingandchemicale

27、ngineering D.,-Alcazar A I, Sharman K. Evolving recurrent neural network architectures by programmingJ.GeneticProgramming,1997:89- 10 Pandya A S, Kulkarni D R, Parikh J C. Study of time series prediction under ernationalSocietyforOpticsandPhotonics,1997:116-11 Studer L, Masulli F. Building a neuro-f

28、uzzy system to efficiently forecast chaotic time Nuclear Instruments and Methods in Physics Research Section A: Accelerators, DetectorsandtedEquipment,1997,389(1):264-12 Rosipal R, Koska M, Farkas I. Chaotic time-series prediction using resource-allocating networksJ.13 Stasinopoulos D M, Rigby R A.

29、Generalized additive s for location scale and inRJ. Journal of isticalSoftware,2007,23(7):1-14 Poon S H, Granger C W J. Forecasting volatility in l markets: A reviewJ. Journal erature,2003,41(2):478-15 ，岳红.随机系统输出分布的建模，控制与应用J. 控制工程2003,10(3): 193-16 Yoon B J, functionsJ.Signaln P P. A multirate DSP f

30、or estimation of discrete probability sing,IEEEionson, 2005,53(1): 252-17 Kedem B, Fokianos K. Time Series Following Generalized Linear Time sJ. s 18 Wills, T.B. Schn, L. Ljung, and B. Ninness, Identification of HammersteinWiener Automatica,2013vol.49,no.1,pp.7019 Akaike H. A new look at the istical

31、 on,1974, 19(6): 716-identificationJ. Automatic Control,IEEE 20 SchwarzG.Estimatingtheofa J.The annalsof istics,1978, 6(2): 461-21 Dai L, Zanobetti A, Koutrakis P, et al. tions of Fine Particulate Matter Species with in the United es: A Multicity Time- ysisJ. Environmental 22 Tibshirani R. shrinkage

32、 and selection via the lassoJ. Journal of the Royal Society.SeriesB(Methodological),1996:267-23 Zou H. The adaptive lasso and its oracle propertiesJ. Journal of the American istical 2006,23 Zou H. The adaptive lasso and its oracle propertiesJ. Journal of the American istical 2006,101(476): 1418-24 K

33、onishiS,KitagawaG.InformationcriteriaandisticalingM.Springer,25 Nagelkerke N J D. A note on a general definition of the coefficient of determinationJ. 1991,78(3):691-26 大, 27 Stasinopoulos D M, Rigby B A, Akantziliotou C. gamlss: Generalized Additive s for ScaleandJ.R package, 2009: 2.0-28 HastieTJ,

34、TibshiraniRJ.GeneralizedadditivesM.CRCPress,29 来生.数理统计： 科30 Parker E N. Dynamics of 1958,128: lanetary Gas and Magnetic FieldsJ. The Astrophysical ARX模型(时变data1=read.table(nasa20120901-20130831.txt,header=T) data2 = data1,c(4:7)mARX模型(时变data1=read.table(nasa20120901-20130831.txt,header=T) data2 = da

35、ta1,c(4:7)m=data2,4 hist(m ) a=var(m) b = len1=length(m) bic.s=numeric() con=gamlss.control(trace=FALSE) LinkFunction1=function(col1,num=length(col1) #v3 = vLinkFun=paste(v1,col11,sep=) if(num 1)for(i1inv_temp=paste(+,v1,col1i1,sep=LinkFun =inkFun,v_temp,#col1function=functionst1=end1 = st1+i-1 col1

36、=c(st1:end1) #Mfunction=M=matrix(1,nrow=len2,ncol=lag1) for (i1 intemp1=lag1+1-i1 temp2 = len1-i1v=mc(temp1:temp2) M,i1 = vlag1 = 50temp1=lag1+ y=mc(temp1:len1) len2 = length(y)m = data2,4 M = cbind(y,M)lag1 = 50temp1=lag1+ y=mc(temp1:len1) len2 = length(y)m = data2,4 M = cbind(y,M)Mas.data.frame(M)

37、 for (i2 inst1=end1 = st1+i2-1 col1=c(st1:end1)Link_mean = LinkFunction1(col1, V) Link_mean=paste(y,Link_mean,sep=)res_ARX_temp = gamlss(eval(parse(text=Link_mean), sigma.fo=1, data = M, family = NO, = con)#bic.ar=res_ARX_temp$sbc bic.s i2 = bic.arin(bic.ssummary(bic.ss=whi i2 = sin(bic.s col1 = col

38、1function(lag1,s,i2) Link_mean=LinkFunction1(col1,V)Link_mean=res_ARX_temp = gamlss(eval(parse(text=Link_mean), sigma.fo=1, data = M, family = NO, control = con)# familyfor(i3in st1=end1 = st1+i3-1 col1=c(st1:end1)Link_sigma=LinkFunction1(col1,Link_sigma=paste(sigma.fo=,Link_sigma,res_ARX_temp = gam

39、lss(eval(parse(text=Link_mean), eval(parse(text=Link_sigma), data = M, family = NO, control = con)# familybic.ar=res_ARX_temp$sbc bic.s i3 = bic.ars= i3 = col1 = col1function(lag1,s,i3) Link_sigma=LinkFunction1(col1,V)Link_sigma=paste(sigma.fo=,Link_sigma,i3 = col1 = col1function(lag1,s,i3) Link_sig

40、ma=LinkFunction1(col1,V)Link_sigma=paste(sigma.fo=,Link_sigma,res_ARX_temp = gamlss(eval(parse(text=Link_mean), eval(parse(text=Link_sigma), data = M, family NO,control=con)#family e = residuals(res_ARX_temp) y_hat=predict(res_ARX_temp) for(l in # y = y_hat=as.ts(y_hat) e = as.ts(e)ts.plot(y,y_hat,e

41、,gpars=list(xlab=S # t 分布les,ylab=DstValue,lty=c(1:3),col=c(1:3)data1=read.table(nasa20120901-20130831.txt,header=T) data2 = data1,c(4:7)m=data2,4 hist(m ) a=var(m) b = len1=length(m) bic.s=numeric() con=gamlss.control(trace=FALSE) LinkFunction1=function(col1,num=length(col1) #v3 = vLinkFun=paste(v1

42、,col11,sep=) if(num 1)for(i1inv_temp=paste(+,v1,col1i1,sep=LinkFun =inkFun,v_temp,#col1function=functionst1=end1=#col1function=functionst1=end1=st1+p-col1=c(st1:end1) #Mfunction=M=matrix(1,nrow=len2,ncol=lag1) for (i1 intemp1=lag1+1-i1 temp2 = len1-i1v=mc(temp1:temp2) M,i1 = vlag1 = 50temp1=lag1+ y=

43、mc(temp1:len1) len2 = length(y)m = data2,4 M = cbind(y,M)M = #for (i2 inst1=end1=st1+i2-col1=Link_mean = LinkFunction1(col1, V) Link_mean=paste(y,Link_mean,sep=)res_ARX_temp=gamlss(eval(parse(text=Link_mean),sigma.fo=1,nu.fo=1,data=M,family=TF, control = con)# familybic.ar=res_ARX_temp$sbc bic.s i2 =bic.ars=whi i2 = scol1 = col1function(lag1,s,i2) Link_mean=LinkFunction1(col1,V)Link_mea