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,?；谝陨戏治龆x() () () = =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