文稿系統(tǒng)辨識(shí)

1、大連理工大學(xué)系統(tǒng)辨識(shí)課題目名稱：Generalized additive sforlocationand(GAMLSSin評(píng)學(xué)大連理工大學(xué)系統(tǒng)辨識(shí)課題目名稱：Generalized additive sforlocationand(GAMLSSin評(píng)學(xué)院（系電專業(yè)儀器儀表學(xué)位類別學(xué)術(shù)型 專業(yè)型 姓名號(hào)： 導(dǎo)師秦完成日期: 2015年6月81 研究背景及意時(shí)間序列(timeseries)1.1 1 研究背景及意時(shí)間序列(timeseries)1.1 圖地磁擾動(dòng)強(qiáng)度Dst指數(shù)時(shí)間序時(shí)間序列分析(Time1.2 描建決策和控1.2 1.2 描述：用圖形(1.1)來(lái)描述時(shí)間序列，或應(yīng)用樣本的自協(xié)方差、

2、樣本相關(guān)即是指從Dst陽(yáng)風(fēng)速度和等離子密度等變量呈現(xiàn)了整數(shù)型泊松分布或?qū)崝?shù)型系統(tǒng)中， 因此，開(kāi)展非平穩(wěn)時(shí)間序列的研究具有重要的理論意義和應(yīng)用價(jià)值產(chǎn)生時(shí)間序列的系統(tǒng)看作隨機(jī)系統(tǒng)，并將所獲得的數(shù)據(jù)集 =| = 0,1,2,的觀測(cè)值量1.2 描述：用圖形(1.1)來(lái)描述時(shí)間序列，或應(yīng)用樣本的自協(xié)方差、樣本相關(guān)即是指從Dst陽(yáng)風(fēng)速度和等離子密度等變量呈現(xiàn)了整數(shù)型泊松分布或?qū)崝?shù)型系統(tǒng)中， 因此，開(kāi)展非平穩(wěn)時(shí)間序列的研究具有重要的理論意義和應(yīng)用價(jià)值產(chǎn)生時(shí)間序列的系統(tǒng)看作隨機(jī)系統(tǒng)，并將所獲得的數(shù)據(jù)集 =| = 0,1,2,的觀測(cè)值量就是要根據(jù)觀測(cè)到的數(shù)據(jù)集建立關(guān)于的條件概率密度模型(|AI,) 。其中，AI

3、于數(shù)據(jù)集國(guó)內(nèi)外研究及發(fā)展現(xiàn)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圖通過(guò)圖 1.3 AR ARMA(AutoRegressive and Moving Average ARARX(AutoRegressive exogenous) ARMA 通過(guò)圖 1.3 AR ARMA(AutoRegressive and Moving Average ARARX(AutoRegressive exogenous) ARMA ARXMA(AutoRegressiveexogenousandMovingAverage

5、況下，圖1.3中的1.4 當(dāng)狀態(tài)變量不可觀測(cè)時(shí)通常用圖1.4所示的狀態(tài)空間方程來(lái)描述時(shí)間序列中，為獨(dú)立同分布的。雖然如HammersteinWiener模型18AR 類模型中要求數(shù)據(jù)的方差為穩(wěn)態(tài)的這一問(wèn)題，Robert Engle 歸條件異方差(AutoRegressiveConditionalHeteroskedasticity2 = 0 + 其中有效的解決一些實(shí)際問(wèn)題。ARCH 2003 經(jīng)濟(jì)學(xué)獎(jiǎng)的計(jì)量經(jīng)濟(jì)學(xué)成果之一。但由于其規(guī)定了此在理論與實(shí)際應(yīng)用上仍具有一定的局限性262 = 0 + 其中有效的解決一些實(shí)際問(wèn)題。ARCH 2003 經(jīng)濟(jì)學(xué)獎(jiǎng)的計(jì)量經(jīng)濟(jì)學(xué)成果之一。但由于其規(guī)定了此在理論與

6、實(shí)際應(yīng)用上仍具有一定的局限性26近任意的分布15。在信號(hào)處理領(lǐng)域，Yoonn 出了類似的方案16。由于核密度估計(jì)模型充分考慮了隨機(jī)輸出變量概率分布的非對(duì)稱 較大。為解決這一問(wèn)題，本研究將基于廣義線性模型(Generalized linear ，GLM模型17 。在統(tǒng)計(jì)學(xué)上，GLM 模型是一種受到廣泛應(yīng)用的線性回歸模式。GLM 模型是基于指 logistic 回歸以及泊松回歸等經(jīng)典的回歸模型。在某些條件下，GLM 究非平穩(wěn)的時(shí)間序列數(shù)據(jù)進(jìn)行建模。AIC(Akaike information criterion)19Information Criterions)201.2所示。此后，Konishi

7、GIC24，GIC 外，Lasso Zou 等還將1-及2-Elastic net 當(dāng)?shù)恼齽t化參數(shù)，Lasso 度以及信息量準(zhǔn)則問(wèn)題也受到一定的關(guān)注。除了1-外，Lasso Zou 等還將1-及2-Elastic net 當(dāng)?shù)恼齽t化參數(shù)，Lasso 度以及信息量準(zhǔn)則問(wèn)題也受到一定的關(guān)注。除了1-等1/2-Bridge 表AIC 、BIC對(duì)(Goodness of Fit)的一本文主要內(nèi)容及創(chuàng)圖1.5 計(jì)出的模型利用決定系數(shù)進(jìn)行評(píng)價(jià)，并利用BC信息量準(zhǔn)則進(jìn)行估計(jì)，選用BIC的原因?qū)⒃诘谌陆o出；最后，將研究方法用于空間天氣預(yù)報(bào)，并對(duì)所提方法進(jìn)行反饋修正。圖1.5 計(jì)出的模型利用決定系數(shù)進(jìn)行評(píng)價(jià)，

8、并利用BC信息量準(zhǔn)則進(jìn)行估計(jì)，選用BIC的原因?qū)⒃诘谌陆o出；最后，將研究方法用于空間天氣預(yù)報(bào)，并對(duì)所提方法進(jìn)行反饋修正。本研究的主要?jiǎng)?chuàng)新在于將基于GLM所模型可以對(duì)決定概率分布性質(zhì)的數(shù)學(xué)期望，方差，偏度(Skewness)與峰度分布，而根據(jù)本研究的結(jié)果發(fā)現(xiàn)，使用基于2 回歸意義下的GAMLSS 模GLM 時(shí)的情況，對(duì)于像服從于分布的時(shí)間序列，有數(shù)學(xué)期望，方差度三GLM GLM 、Akantziliotou 、Stasinopoulos 2002 年提 additive s)28的2 回歸意義下的GAMLSS 模GLM 時(shí)的情況，對(duì)于像服從于分布的時(shí)間序列，有數(shù)學(xué)期望，方差度三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 經(jīng)過(guò)研究 發(fā)現(xiàn)，GAMLSS是一種(半獨(dú)立式)參變的回歸模型，這意味著假設(shè)響GAMLSS 中對(duì)于響應(yīng)變量要服從指數(shù)族分布這一要求是很變量或是隨機(jī)影響都是可以的。因此，GAMLSS 特別適合那些 GAMLSS模型中，如果令 = , = ( , , , ) ， 1211)，1是= ( )的逆如果 是單

10、數(shù)的話， 將服從與exp( 成 2比例的錯(cuò)誤的先驗(yàn)密度函數(shù)。由公式(2.1)GAMLSS GAMLSS () ( ) = = + ( 其中 = 1,2,3,4時(shí)分別為, , , ，()為關(guān)于解釋變量未知函數(shù)，為在GAMLSS()1() = = 如下式所示，模型(2.5)可以擴(kuò)展為具有 ( ) = = ( , ) + 如下式所示，模型(2.5)可以擴(kuò)展為具有 ( ) = = ( , ) + ( (3.26)稱為非線性半?yún)?shù)累加模型。如果 = 可以得到非線性參數(shù)GAMLSS() = = (,如果() = (2.7)與(2.5)GAMLSS (3.27)中的() GAMLSS 模型可以看成線性和非線

11、性參數(shù)模型的結(jié)合。因此 稱(2.7)或(2.5)GAMLSS 模型。在GAMLSS 模型中，本研究利用罰似然函數(shù)(2.8)來(lái)進(jìn)行模型參數(shù)與 = 1 2其中=log (|)。在本研究GAMLSS3 進(jìn)一步研究設(shè)3.1 經(jīng)典時(shí)間序列模表3.13.1 AR可以噪聲較為復(fù)雜的ARMA在時(shí)間序列分析中很多模型都可以的用狀態(tài)空交替計(jì)算 = 1(1 )和 = (3 進(jìn)一步研究設(shè)3.1 經(jīng)典時(shí)間序列模表3.13.1 AR可以噪聲較為復(fù)雜的ARMA在時(shí)間序列分析中很多模型都可以的用狀態(tài)空交替計(jì)算 = 1(1 )和 = (濾波多時(shí)間序列分析都可以轉(zhuǎn)換化為狀態(tài)空間模3.2 GLM模獨(dú)立。在數(shù)據(jù)處理過(guò)程中認(rèn)為響應(yīng)的過(guò)

12、去值的條件期望 = E|1是一個(gè)關(guān)于協(xié)變?yōu)轫憫?yīng)的過(guò)去值的條件期望 于躍遷概率為 的二元時(shí)間序列來(lái)說(shuō)， 對(duì)協(xié)變量進(jìn)行線性回歸分析 = 的Poisson或者Gamma 模型6(3.1)，這是一種嚴(yán)格服從系統(tǒng)性和隨機(jī)1() = 2() = 這里 = 1, ( ( ; ,) = exp+ ( ;這里 = 1, ( ( ; ,) = exp+ ( ;其中 () ， 離散參數(shù)是已知的權(quán)值或者過(guò)去的權(quán)值， 是分布的自然2 系統(tǒng)性部分：對(duì)于 = 1,有單調(diào)函數(shù)g(.() = = = 函數(shù)g(. )是當(dāng)變量向量 接函數(shù)。當(dāng)一維響應(yīng)向量 = 0 + 11 + 22 + 3cos 或 = 0 + 11 +22 +

13、31 + = + )+ ) 其中, 1 = ,1(1),(),1(1), = (,1,1,() = = + ) ) + 其中= () ，因?yàn)榘俗曰貧w和滑動(dòng)平均兩個(gè)部分，被稱為GLARMA(Generalized Linear Autoregressive Moving Average)GARMA(Generalized Autoregressive Moving Average)。() = = 目前Box-Cox = 1 =0 表示對(duì)數(shù)線性模型。需要注意的是當(dāng)未知參數(shù)為相乘關(guān)系時(shí)，在估計(jì)GARMA模型時(shí)在上述分析中，關(guān)于1， = E|1 的解釋如下。假設(shè)= ( )序列() = = 目前Bo

14、x-Cox = 1 =0 表示對(duì)數(shù)線性模型。需要注意的是當(dāng)未知參數(shù)為相乘關(guān)系時(shí)，在估計(jì)GARMA模型時(shí)在上述分析中，關(guān)于1， = E|1 的解釋如下。假設(shè)= ( )序列p 維向量 ( = 1,)。由1，1確定的11 = (1,2 ,1,2 知道對(duì)于 = 1,有(|1 = 1同時(shí)對(duì)： ( (; ,) = = 將等式(3.8)兩邊同時(shí)對(duì)：(;,|1) (;,|1) (; , |1) = = ( ( (;,|1) = (;,|1)(; , |1) = () (; , |1) =：Var|1 = 因?yàn)閂ar|1 0，是單調(diào)的，因此可以將(3.10) = 因此由連接函數(shù)(因?yàn)閂ar|1 0，是單調(diào)的，因

15、此可以將(3.10) = 因此由連接函數(shù)( ，1 ()1可以導(dǎo)出 = = 通過(guò)以上分析它的特點(diǎn)如表3.1表在統(tǒng)計(jì)學(xué)上，GLM此模型假設(shè)實(shí)驗(yàn)者所量GLM測(cè)的隨量的分布函GLM數(shù)與實(shí)驗(yàn)中系統(tǒng)性效應(yīng) (即非隨機(jī)的效應(yīng))可由一 連 接 函 數(shù) (link function)建立起可以解GLMGLM 模型服從條件概率密度函數(shù)為(|21AI)的概率分布。其中AI為輔助信息量，然而在大多數(shù)情況下 AI AI ) = = (|1,2,1)(,1,2,= (|1,2,1)(1|1,2,2)(,1,2,在有協(xié)變量參與時(shí)似然函數(shù)() = (1122) 其中 = (112 其中 = (112211)， = (1122

16、11)() = (; , ()=log(;, =+( ; ( ) ( = +(;其中() ()1 = 11()，并且從(3.9)可以推出 = ( ,1 = 其中 = 1,。由(3.9)和(3.10) = 1( Var 又因?yàn)? = = Var|1其中 = 1,?；谝陨戏治龆x() = = Var|1其中 = 1,?；谝陨戏治龆x() () () = =1 其中2() = Var|1。對(duì)于 = 1,定義() () = =1 2() E = 1 2() 假設(shè)E() = 0，當(dāng) () () = E() = log() = 由此可以解出AR 模型中，通常定義其噪聲ARGAMLSSDst數(shù)據(jù)進(jìn)行模型

17、估計(jì)與選擇時(shí)，首先假設(shè)其服分布利用AR 模型進(jìn)行數(shù)據(jù)的建模3.1 為方差2分布概率密度函數(shù)曲線。如果(量1() ( 由此可以看分布是典型的指數(shù)族分布。由圖3.1 可以看曲線關(guān)于軸對(duì)稱，并且|3.1 為方差2分布概率密度函數(shù)曲線。如果(量1() ( 由此可以看分布是典型的指數(shù)族分布。由圖3.1 可以看曲線關(guān)于軸對(duì)稱，并且|的值越大，()的值越趨近于零，當(dāng)| 時(shí)() 0 AR 模型中定義噪聲(0,2)望和方差，因此在進(jìn)行模型的估計(jì)與選擇是只需估計(jì)和。GAMLSS = +1 log()=+0基于AR由于GAMLSS服從分布時(shí)進(jìn)行對(duì)數(shù)學(xué)期望，標(biāo)準(zhǔn)差度稱設(shè)(0,1)，2且和 度為 的 分布29。 其中，

18、 如果是 變量， 其分布稱為12 i.i.d.(0,1) =度為的2變量，其分布稱度為的2分布，記為2 度為 的 分布29。 其中， 如果是 變量， 其分布稱為12 i.i.d.(0,1) =度為的2變量，其分布稱度為的2分布，記為2量，其概率密度函數(shù)可以寫(xiě)為(3.21)3.2 2) 2(1n3.2 3.1 分布曲線很相似。二者均為關(guān)于原點(diǎn)對(duì)稱，單峰值的偶函數(shù)，并在 = 0但分布相比需要注意的是，()當(dāng)且僅當(dāng) 1) + )22n1() ( )( 當(dāng) 2時(shí)， + )22n1() ( )( 當(dāng) 2時(shí)，() = 0；當(dāng) 3時(shí)，() 當(dāng) = 1時(shí)，t11() = (1 + 其中 ()當(dāng) 時(shí)，變量的極限

19、分布為(0,1)Dst GAMLSS = +1 log( ) = +2 0 log() = +3 0模型的評(píng)價(jià)與選determination)(3.22) 2 = 1 決定系數(shù)2是相關(guān)系數(shù)的平方， 數(shù)據(jù)對(duì)隨機(jī)模型的擬合程度的一個(gè)隨機(jī)變量。決定系數(shù)通常基于相關(guān)信息對(duì)輸出值進(jìn)行 或者對(duì) 進(jìn)行檢驗(yàn)。目前根據(jù)不在以上兩種情況中，2 (0,1)。當(dāng)性回歸模型，為解決這一問(wèn)題，Cox 和 Snell 在，()上節(jié)給出了當(dāng)自回歸階數(shù)給定情況下模型參數(shù)數(shù)AIC BIC 來(lái)確定最優(yōu)的自回歸階數(shù)(延時(shí)) = 1, BICBIC 最小時(shí)的4由上一章的內(nèi)不難發(fā)現(xiàn)，本研究是將GLM 模型擴(kuò)展到時(shí)間序列對(duì)非 風(fēng)與Dst指

20、數(shù)的相上節(jié)給出了當(dāng)自回歸階數(shù)給定情況下模型參數(shù)數(shù)AIC BIC 來(lái)確定最優(yōu)的自回歸階數(shù)(延時(shí)) = 1, BICBIC 最小時(shí)的4由上一章的內(nèi)不難發(fā)現(xiàn)，本研究是將GLM 模型擴(kuò)展到時(shí)間序列對(duì)非 風(fēng)與Dst指數(shù)的相風(fēng)與Dst指數(shù)介Dst(Disturbance Storm Time)Worldenter eomagnetism, Kyoto 2012的4.1 的左側(cè)，右側(cè)為4.1 Nagelkerke 在1991 年對(duì)其有了進(jìn)一步的研究25，由于關(guān)于非線性模型22 = 1 ()11 風(fēng)的速度一般為200800km/s極風(fēng)一詞是由尤金派克301950 1960 于等、4.2Dst4.2 Dst D

21、st工Dst仿真實(shí)驗(yàn)結(jié)11 風(fēng)的速度一般為200800km/s極風(fēng)一詞是由尤金派克301950 1960 于等、4.2Dst4.2 Dst Dst工Dst仿真實(shí)驗(yàn)結(jié)4.3DstDstDst指數(shù)數(shù)據(jù)的曲4.3 4.4 為數(shù)據(jù)的直方圖。4.4Dst4.4Dst由于 分布的似然函數(shù)較為簡(jiǎn)單，并且容易計(jì)算。因此本研究首先考慮數(shù)據(jù)服從4.4 100 的分析了解到分布有期望、方差、 度這三個(gè)對(duì)數(shù)據(jù)分布情況影響較大的統(tǒng)計(jì)量。在本研究中嘗試對(duì)期望、方差、 度分別進(jìn)行動(dòng)態(tài)估計(jì)，其具體過(guò)程如下：假設(shè)數(shù)據(jù)的方差度為定常，對(duì)期望進(jìn)行動(dòng)態(tài)建模假設(shè)數(shù)據(jù)度為定常，分別對(duì)期望、方差進(jìn)行動(dòng)態(tài)建模假設(shè)數(shù)據(jù)的期望、方差度均為時(shí)變的

22、，分別對(duì)三者進(jìn)行動(dòng)態(tài)建模本研究中利用R中的GAMLSS數(shù)據(jù)包計(jì)算地磁Dst指數(shù)時(shí)間序列的模型參數(shù)、決定系數(shù)以及BIC 選擇結(jié)果。方差為定常分布數(shù)據(jù)的仿真實(shí)本仿真實(shí)驗(yàn)中，令最大延時(shí) = 504.5AR4.5Estimate 4.5AR4.5Estimate 為模型參數(shù)的估計(jì)值，Std.Error SBCBICBIC值。通過(guò)仿真結(jié)果可以看出，BIC選擇出的最優(yōu)的最大延時(shí) = 34.6 為模型擬合后的曲線。圖AR模型擬合曲4.7 圖時(shí)變AR模型的仿真結(jié)BIC 信息量準(zhǔn)則選擇的最優(yōu)的最大延時(shí) = 圖時(shí)變AR模型的仿真結(jié)BIC 信息量準(zhǔn)則選擇的最優(yōu)的最大延時(shí) = 4.2.1 BIC 4.8 圖時(shí)變的AR

23、模型擬合曲方差為常數(shù)的服從于4.9 4.9 分布仿真結(jié)果(方差定常于4.9 4.9 分布仿真結(jié)果(方差定常 = 504.9度都是常數(shù)的情況下，BIC4.10 4.10 分布擬合曲線(方差定常度為定常的服從。由4.2.3 節(jié)。由4.2.3 節(jié)度都是常數(shù)，通過(guò)BIC 選擇出的圖4.11 為估計(jì)結(jié)果。圖4.11 分布的仿真結(jié)果度定常BIC對(duì)數(shù)據(jù)的方差進(jìn)行估計(jì)時(shí)，BIC 信息量準(zhǔn)則的選擇結(jié)果 = 2，可以認(rèn)為這是一個(gè)較為4.12 為擬合曲線。4.12 服從在及 = 504.13 時(shí)變的在及 = 504.13 時(shí)變的4.13 BIC選擇出的最優(yōu)延時(shí)為3BIC擇結(jié)果為 = 2BIC 選擇出的最優(yōu)延時(shí)為 =

24、 24.14 仿真小BIC4.1 3.4.1 BICAR33AR32仿真小BIC4.1 3.4.1 BICAR33AR323224.1 AR 對(duì)于服從分布的模型，只有期望時(shí)變的這種情況，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,安，杜金觀，潘一民時(shí)間序列的分析與應(yīng)用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ī)系統(tǒng)輸出分布的建模，控制與應(yīng)用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 來(lái)生.數(shù)理統(tǒng)計(jì)： 科30 Parker E N. Dynamics of 1958,128: lanetary Gas and Magnetic FieldsJ. The Astrophysical ARX模型(時(shí)變data1=read.table(nasa20120901-20130831.txt,header=T) data2 = data1,c(4:7)mARX模型(時(shí)變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