誕生這篇論文的緣起與波折
當初進來北大時，我最想要找的指導教授就是林財川老師，因為老師的研究專長是時間序列，是我想專攻的領域。當時跟老師面談時，老師跟我說的研究主題其實比較偏向 Hammerstein system、Projection Pursuit Regression(PPR)跟neural networks (NN) 跟一般 time series model 的預測能力比較，當時聽得似懂非懂。Hammerstein system 不是一般時間序列課本會出現的東西，而且研究它的學者很喜歡把直接模型結構畫成流程圖，對當時的我來說很魔幻。1老師當時還說像 Hammerstein system、Projection Pursuit Regression(PPR)、neural networks (NN) 都是非參數方法(Nonparametric Methods)，當時也不知道這是甚麼。2
因為當時的研究偏向預測的關係，老師讓我去找了真實世界多維度的資料來做 case study，於是我考慮了幾個方向後，決定使用沒有人做過的台中捷運每日流量、測站溫度跟濕度，並且與老師討論後，使用 ARIMA 跟 SARIMA 作為一般 time series model 的對照組。但過了一段時間，到了碩一下的時候，老師突然跟我說方向要改，因為隔壁有位老師指導學生也用了類似的題目，也就是說我們撞題目了！
雖然撞題目，不過老師想過以後說，要使用的模型完全不用改，只要改成從有效自由度(effective degrees of freedom)的角度出發。這個角度比純粹的預測還要特別，讓我覺得我的題目跟別人與眾不同！看了非常喜歡，之後大方向就定案了。不過也因為角度變了，中捷的資料也就沒用了。
至於非參數方法，後來我發現非參數的定義是模型沒有來自機率分布，但這個定義在實務上不夠明確，可能會發生 A 學者認為這個模型是非參數方法、B 學者卻不這麼認為的情況。3於是後來我改成線性跟非線性模型的比較，只要模型有非線性的部份就會被定義成非線性模型。
另外，在過程中，原先預定採用的 Hammerstein system 先是被替換成 Wiener system，兩個東西系出同源，只是結構稍有不同。後來因為這個流派的方法不會做估計這段，所以後來沒有採用。我看了這些方法的 paper 將近一年，太可惜了 QQ，有機會再來講這是甚麼東東。
有效自由度是甚麼
歷史背景
要說明有效自由度的概念，就需要先了解自由度這個東西。自由度的概念在商統中就有學過，但對初學者來說仍是需要花費一點時間才能真正理解的東西。這也難怪，自由度的概念最初只是存在天才腦海中的東西，文獻上只有一句話帶過，沒有正式定義。4有一段時間要怎麼向學生正確解釋自由度這個概念還成為學者之間的討論話題，(Walker, 1940)就有提到要如何正確解釋自由度，他的版本近似於大學修迴歸分析時的老師交給我的版本。5自由度也可以再進一步歸類，如果是代表整個模型的自由度，將模型寫成矩陣形式，在線性迴歸模型中自由度的概念會像這樣：
Consider a linear model \[\mathbf{Y}=\beta \mathbf{X} +\varepsilon\] where \(X\) is a \(n \times p\) matrix, \(\beta\) is a \(p \times 1\) matrix, and \(Y\) is a \(n \times 1\) matrix. \[\hat \beta = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}\] and the fitted value \(\mathbf{\hat Y}=(\hat y_1,\dots \hat y_n)\) is \[\mathbf{\hat Y} =\mathbf{X}\hat\beta = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}\]
此時自由度 \(p\) 就是來自以下 hat matrix(\(\mathbf{H}\))的跡(trace)：
\[\mathbf{H}=\mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\] \[\mathbf{\hat Y} =\mathbf{H}\mathbf{Y}\] \[tr(\mathbf{H})=p\]
上面提到的 \(\mathbf{\hat Y} =\mathbf{H}\mathbf{Y}\) 叫做 linear smoother ，當估計模型的估計值跟真實觀測值可以直接寫成這樣時，我們可以說兩者的關係線性可微。可想而知，能夠直接這樣寫的都是線性模型。
但，自由度的概念因為太抽象，加上在線性模型中它的數值等價於模型使用的參數數量，所以一般在應用上只會用參數數量來簡單表示。
可是，這樣真的好嗎？
只用「參數數量」表示自由度會有的問題
舉個例子，今天我使用今日的雨量(\(X_{t}\))跟前一天的雨量(\(X_{t-1}\))來預測當日捷運站的日進入人次(\(Y_{t}\))。我建了簡單的模型 1 號：
\[Y_{t}=\beta_1X_{t}+\beta_2X_{t-1}+\varepsilon\]
跟稍微有點變化的模型 2 號：
\[Y_{t}=\beta_1(X_{t})^{0.5}+\beta_2X_{t-1}+\varepsilon\]
試問，這兩個模型的自由度是否相同？
如果只用參數數量表示，那兩者自由度是相同的，但是模型 2 號卻比模型 1 號多個對 \(X_t\)開根號的步驟，直覺上來說很奇怪。如果回歸到自由度的數學定義，那模型 2 號會因為沒有辦法寫成 \(\mathbf{\hat Y} =\mathbf{H}\mathbf{Y}\) 而沒有定義，也就是說，模型 2 號的自由度無法計算！
這時候就需要我們有效自由度了
有效自由度的概念最初是由(Ye, 1998)提出，背後隱含模型複雜度的概念，也就是「實際上限制模型發展的參數量」，這明顯比參數數量更能反映自由度的理念。它的概念是從自由度的概念延伸，可視作一般化版本。它的基本概念就是對現有模型的估計數值作擾動，只要這個估計模型穩健，那對擾動後的數值與原有的估計做泰勒展開，就可以成功拼湊出 \(\mathbf{H}\)來
\[(\widehat{\mathbf{Y} + \Delta}) - \hat{\mathbf{Y}} \cong \mathbf{H} \Delta\]
然後再對\(\mathbf{H}\)取 trace 即可，聽起來非常簡單。
可是…
它的解實際上並不好找，因為很多模型的\(\mathbf{H}\)是得不出來的，必須要直接求 hat matrix 的對角元素 \(h_{ii}\)再相加，在這過程中，會涉及到三個角色：由真正模型( true model)產生的自然資料(observation value)、最適合估計它的估計模型產生的估計值(true model estimation)跟現有的估計模型產生的估計值(model estimation)。想當然，我們根本不知道最適合估計它的估計模型的估計會長甚麼樣，GG。
那就用數值方法來模擬吧
為了解決這個問題，(Ye, 1998)就用蒙地卡羅的方式來估計有效自由度。要達成這個目的，需要設定最一開始擾動參數的變數程度(用 \(\tau\) 來表示)。在文獻中，它主要應用在 Tree 類的迴歸模型上，並初步得出當模型的變異數是\(\sigma^2\)，\(\tau\)的設定範圍在 \([0.5\sigma,1\sigma]\) 時估計結果相對穩定，\(\approx 0.6\sigma\)則是合理的選擇。
有效自由度的應用與疑問
自由度雖然是個不起眼的小東西，但是現在常用的模型選擇標準(model selection criteria)卻會直接或是間接用到它，例如 AIC 和 Cross-Validation(CV)等，這些方法現今常用在更為複雜的模型上，從中選出表現最好的模型。如果繼續使用簡單的參數數量來表示自由度，很可能會連帶削弱模型選擇標準的效果，進而導致不小心選出泛化能力不佳的模型，簡單的說，能好好預測未來的模型反而會因此被埋沒。有效自由度便能克服這樣的問題，而當它用在 Cross-Validation 上就叫做 Generalized Cross-Validation (GCV)。
不過，有關有效自由度的直接應用反而不多，如前所述，這是一個有點難直接求解的東西，也涉及到最一開始擾動參數的變異程度設定。就我所知，目前除了 Lasso 跟 B-Spline 等模型有快速的公式解可以直接得知有效自由度外，一般模型的有效自由度皆需要用模擬的方式得出。但是最初在(Ye, 1998)的文獻中只有用在 Tree 模型上，他所謂合理的參數設定是否適合用在時間序列模型上，是值得驗證的問題。
本次用到的模型
相較於有效自由度，我在這裡用到的模型反而滿普通的 XD，只有 ARIMA、SARIMA、Projection Pursuit Regression (PPR)跟 Neural Networks (NN)。其中 NN 是以前神經網路還沒有衍生出很多模型時使用的稱呼，現今在機器學習的領域裡通常指的是前饋式神經網路(Feedforward Neural Network, FNN)。6因為在時間序列領域仍有學者稱呼 NN ，我的老師也這樣說，所以我也沿用這個稱呼。
PPR 也是個有趣的東西，它最早是由 Friedman 在 1981 年提出(Friedman & Stuetzle, 1981)，一度成為能跟 NN 交鋒的對象。如果有讀過 Friedman 的 paper，會發現他的研究主要關注讓模型自動找出最佳估計， PPR 就是他的研究成果之一。它的模型結構設定的很 general ，效仿 Projection Pursuit 的精神自動找出最佳投影。它的結構與 NN 有幾分相似之處，但是在後來學者討論指出這個模型的泛化能並不好，本次我的研究也可以藉由有效自由度清楚發現這一點。
至於線性模型的代表 ARIMA 跟 SARIMA，我就懶得說了，因為時間序列都會介紹到這兩位，就不特別說明了。
研究成果與感受
我的研究架構最初是沒有跑模擬，但是因為前面產生的疑問，讓我想跑模擬驗證，老師也准了。結果…
跑模擬真的好痛苦!!!
我已經有點忘記跑模擬到底卡了多久，但我深深的記得那個過程有多痛苦。即使我已經設定最簡單的模型，但是關於實驗該如何設計卻是我從來沒有想過的問題，也不知道要如何以最好的方式寫下觀察成果。因為這樣我中間修改了很多次，一度看不見論文完成的盡頭。
此外，雖然模擬有效自由度的方法在 1998 年就有了，但是從來沒有人把它包成一個獨立的 package ，因為它涉及不同模型的演算法，除了前面提過可以快速求出的 Lasso 跟 B-Spline 以外，只能自己寫!!雖然老師給我用 S-Plus 寫的 code ，但說實話，跟 R 的邏輯還是有點不同，我沒有辦法直接轉換…
當然我也求助過各種 AI 服務，但結論是他們吐出來的東西邏輯跟文獻似乎有差，我還是不懂也不放心。最後，我決定先照文獻課算法，不會再一小步一小步問 AI 怎麼寫、為什麼要這樣寫 ，總算是成功刻出來了，每個步驟都照著文獻走，我也完全理解它們的涵義。
不過也因為前面的模擬跑了千千萬萬次，我在跑實際資料時花費的時間反而不多，唯一有差的是資料前處理跟處理時間長很多。我這次的 Real data case 改成現在依舊很熱門的太陽黑子數，這是一種有週期性、但不好預測的資料。因為資料有右偏斜的問題，直接建模會出事，所以好的資料轉換是很重要的。雖然我試著照著文獻做資料轉換，不過最後發現表現最好的還是 Tang 在專書(Tong, 1993)中用的 root transformation。7
然後我在使用文獻推薦的模型算有效自由度時，發現有的模型加入擾動後，反而會因此求不出解！這可能是因為時間序列模型也有自己的收斂方法，加入擾動反而會干擾它求解。目前可行解決方式只有把\(\tau\)設小，使擾動的變異度變小，但設太小在模擬結果中顯示算出來的結果會低估，所以這樣也不太好。但我也發現不管是怎麼樣的設定，PPR 的有效自由度很明顯都會高過於所有模型，代表它的模型結構太複雜，影響它的泛化能力；反之，與它結構略顯相似的 NN 的模型複雜度反而低很多，展現出良好的泛化能力。
我做完這個研究後，深刻理解到有效自由度為何很少人在研究，還有它在理論上的美麗之處與現實的骨感。8不過，我也從中發現並且重新審視 Ye 參數選擇的意義，口委也給予我許多寶貴的回饋意見。
總結
我覺得有效自由度絕對是個好東西，它幫助我們理解天才腦中的想法。希望在不久的將來，像我們這樣的凡人也能過算出有效自由度來選出最好的模型。
這裡有出現的文獻
腳註
後來發現 Box 寫的時間序列課本也是把時間序列想成流程圖，也許老師的啟蒙課本是 Box 寫的也說不定。↩︎
比較符合台灣這裡的譯法是意譯的「無母數方法」，表示不用假定機率分布的模型，「非參數方法」反而不夠直覺，不過因為老師都這樣說，所以於予保留。↩︎
舉例來說，像 NN 有人認為它是非參數方法，但有人卻認為它是半參數方法(Semiparametric method，又稱半母數方法)，因為它的模型結構的一些部分仍來自機率分布。↩︎
最初提出這個概念的是費雪，但他甚至沒有為自由度命名，只是用「number of experiments」來帶過。↩︎
我大學的迴歸分析是魏文翔老師教的，也是位好老師！↩︎
不過我自己覺得神經網路家族模型稱呼一直都亂七八糟，合理懷疑不同領域跟學者之間根本沒有嘗試統一稱呼過。↩︎
但這個想法最後還是有被口委 challenge 。↩︎
我做到後面甚至也忍不住像文藝復興時期的科學家跟數學家那樣開始覺得上帝也許真的存在，因為有效自由度的理論解很美麗，但我們這些普通人只能用醜陋的數值方法算出醜不拉基的解(更糟的是它還不一定是正確的)，來榮耀只有上帝知道的理論解。我不排除這是因為模擬跑得太痛苦才這樣想。↩︎