B-Splines不同節點選擇方法之比較 - 政大學術集成

全文

(1)國立政治大學統計學系碩士學位論文. B-Splines 不同節點選擇方法之比較 The comparison between different 治methods of knots selection. 立. 政. for B-Splines. 大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. 指導教授：黃子銘博士研究生：胡子卿撰. 中華民國一百零六年六月.

(2) 要本文以 B-Spline 的框架研究比較兩種不同的節點估計方法。第一種方法是通過最優化特定的目標函數並結合相對應的選擇標準選擇出最優化的節點組合。第二種方法則基於幾何控制多邊形的特性將內部節點的選擇過程與幾何圖形聯繫起來，省去了最優化的過程。另外，本文採用『節點估計時間』與『誤差平方和』（Mean Squared Error）來評價兩種方法的估計結果。通過分析各種不同模擬數據下兩種方法的表現情況，本文的主要發現是：第一，無論哪種資料，第二種方法在計算速度上都是大幅領先第一種方法。第二，在數據資料較小的情況下，第一種方法中由 Lindstrom 提出的算法並不能很好的配飾模型，最後的估計誤差較大。而在數據資料較多的情形下，誤差與其他方法較為接近。第三，第一種方法中沒有懲罰項的算法在所有驗證. 政治大估計方法提供了很具價值的參考信息。立弧線函數. 節點. 估計. ‧. io. sit. y. Nat. n. al. er. ：. 學. 關. ‧ 國. 過的數據中，其表現是所有方法中最穩定且估計誤差最小的。這些發現為如何選擇恰當的節點. Ch. engchi. I. i n U. v.

(3) Abstract This study compares two different methods of knot selection for B-Spline. The first one chooses the best knots through optimizing specific objective functions and corresponding criterion. Based on some properties of geometric control polygon, the second one connects the knot selection process with geometric figures, which avoids the tedious optimization. On the other hand, we use the time for estimation and the mean squared error to evaluate the performance of these two methods. There are three main findings of this study. The first finding is that the calculation speed of second method is much higher than that of the first one. Secondly, the algorithm proposed by Lindstrom in the first method is not stable and its estimation error. 政治大 proposed by Lindstrom becomes better as the sample size increases. Thirdly, the performance 立 is larger when the sample size is small. On the contrary, the performance of the algorithm. of the algorithm without penalty term in the first method is always better than the second. ‧ 國. Spline. knots. ‧. estimation. io. sit. y. Nat. n. al. er. Key Words:. 學. method.. Ch. engchi. II. i n U. v.

(4) 致謝本文是在黃子銘老師的悉心指導下完成的。大到確定選題，研究方向，探索文章思路，選擇分析方法；小到論文規範格式，文本編輯軟體等都融入了黃老師的心血和教導。在黃老師一年半的教誨中，我深深感受到了她嚴謹求實的治學態度，以人為本的研究理念，和藹可親的長者風範，以及對學生無微不至的關懷。所有的這些都讓我不敢有絲毫懈怠，努力做到最好，並暗暗在心中以黃老師的治學之道，為人之道作為我自己的榜樣。我相信，這些將會使我受益終身。在這裏，我想衷心地向黃老師致以無限的敬意和感謝！即將離開台灣了，我還要向政大商學院的其他老師衷心說一聲謝謝，特別是楊素芬老師、鄭宗記老師、陳麗霞老師、江振東老師、鄭宇庭老師、余清祥老師、洪英超老師、薛慧敏老師. 政治大師恩浩蕩，友情無價。我要感謝所有在政大相識的朋友，是你們陪我度過短暫但卻充滿迴立. 等。謝謝你們教給了我各們專業知識，也教給了我許多為人處世的道理。子卿感激不盡。. ‧ 國. 學. 憶的兩年。特別是統計系的各位小夥伴們，你們的睿智，你們的友善無不讓我不捨。我們一起度過學習和生活中的點點滴滴早已內化成我一生中無比珍貴的迴憶。. 我還要感謝我的家人們，感謝你們對我的“放任”，讓我擁有選擇生活的自由。更感謝你們. ‧. 對我方方面面的支持，子卿能走到現在離不開你們給予我的力量。正是你們的付出，才有我一. sit. y. Nat. 步一步攀登求學高峰的機會。. 逝者如斯夫，不捨晝夜。來台的兩年匆匆而去，甚至未能好好地去感受她的美好便要說再. io. n. al. er. 見，實為遺憾。然而人生不就是在一次次分別中找到屬於我們自己人生的真諦嗎？每個人總會. i n U. v. 有屬於自己的那個終點站，但此時此刻，我要在台北站下車了，感恩與我一起分享沿途風景的. Ch. engchi. 人們，你們的故事，你們的一切都是我勾勒美好記憶的點睛之筆。幸運的是，在我即將開啟另外一段旅程的時刻，我遇見了楊貝妮，雖然前方路途遙遠且多有艱險，但內心依舊充滿盼望與堅定，希望與你協手走下去。革命尚未成功，同志仍需努力。就用我最愛的一首詩來結束我在台北的生活吧。 Two roads diverged in a wood, and I—I took the one less traveled by,And that has made all the difference. 感謝日月乾坤，感謝天地萬物，感謝宇宙星辰。. III.

(5) 目. 錄. 要.................................................................... I. ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. II. 致謝 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III 目錄 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. V. 圖目錄 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VI. 政治大. 表目錄 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VII. 立. 1. 1.1. 弧線函數簡介與研究動機 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1. 1.2. 研究目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2. 1.3. 研究框架 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3. 學. 前言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. ‧. ‧ 國. 1. 文獻探討 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3. 研究方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 5. sit. al. er. 最佳化方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. n. 4. io. 3.1. 4. y. Nat. 2. iv n C . . . .h . .e . . . . . . . . . . . .i . . U ngch . . . . . . . . . . . . . . . . . . . . . . . . .. 6 6. 3.1.1. 基於 GCV 和修正後 GCV 的估計 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 6. 3.1.2. 帶懲罰項的估計. 7 9. 3.2. GeD 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3.3. 固定等距節點 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10. 模擬資料分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.1. 資料產生與模擬流程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11. 4.2. 實證結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4.2.1. 節點選擇比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12. 4.2.2. 估計精度比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17. 結論與建議 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25. IV.

(6) 5.1. 研究結論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25. 5.2. 研究建議 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26. 參考文獻 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27. 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. V. i n U. v.

(7) 圖目錄 4.1. σ = 0.1,n = 40,KG = 3 下節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12. 4.2. σ = 0.1,n = 40,KG = 5 節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12. 4.3. σ = 0.1,n = 40,KG = 7 下節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13. 4.4. σ = 0.1,n = 400,KG = 3 下節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14. 4.5. σ = 0.1,n = 400,KG = 5 下節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14. 4.6. σ = 0.1,n = 400,KG = 7 下節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15. 4.7. σ = 0.01,n = 400,KG = 3 下節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15. 4.8. σ = 0.01,n = 400,KG = 5 下節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16. 4.9. σ = 0.01,n = 400,KG = 7 下節點估計數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16. 立. 政治大. 4.10 異常情況下估計表現 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20. ‧. ‧ 國. 學. 4.11 正常情況下估計表現 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20. n. er. io. sit. y. Nat. al. Ch. engchi. VI. i n U. v.

(8) 表目錄 4.1. k G = 3, σ = 0.1, n = 40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17. 4.2. k G = 5, σ = 0.1, n = 40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18. 4.3. k G = 7, σ = 0.1, n = 40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19. 4.4. k G = 3, σ = 0.1, n = 400 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21. 4.5. k G = 5, σ = 0.1, n = 400 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21. 4.6. k G = 7, σ = 0.1, n = 400 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22. 4.7. k G = 3, σ = 0.01, n = 400 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23. 4.8. k G = 5, σ = 0.01, n = 400 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23. 4.9. k G = 7, σ = 0.01, n = 400 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24. 政治大 4.10 沒有誤差的情形下各方法表現情況，n=400 立. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. VII. i n U. v. 24.

(9) 1 前言 1.1 弧線函數簡介與研究動機弧線函數 (Spline) 是一種特殊的函數，在不同的定義域區間上由多項式分段定義。早期 Splines 函數多用於工程製圖方面，後隨著計算機的普及化，其在工業設計方面的應用也越來越豐富。其數學定義可見 [21]，具體如下：假設函數 S 為從區間 [ a, b] 到實數 R 上的映射，即： [ ] S : a, b → R 其次，spline 是分段多項式，其內部分段點稱之為節點，在上述映射的基礎上給定 k 個節點，. 政治大. 這 k 個節點位於區間 [ a, b] 之間迴並有如下關係：. 立 a = γ < γ < ··· < γ < γ 0. k. 1. k +1. =b. ‧. ‧ 國. sit. Nat. Pi : [γi−1 , γi ] → R. y. Pi ，即如下. 學. [ ] 假設函數 S 在不同的區間 γi−1 , γi , 1 ≤ i ≤ k + 1 由於不同的限制條件所產生的多項式函數為. io. er. 則弧線函數具有如下形式：. n. a Sl (γ) =P (γ), γ ≤ γ < γ i v S(γC ) = P ( γ ), γ ≤ γ < γ h. e n g c h i U n .. 1. 0. 1. 2. 1. 2. (1.1). S(γ) = Pk+1 (γ), γk ≤ γ < γk+1 令函數 S 裡面多項式最高次的次數為 m − 1，則 m 稱為該弧線函數的階數。弧線函數還具有在內部節點 γi 處 m − 2 階連續可微的性質，即 ( j). ( j). Pi (γi ) = Pi+1 (γi ), j = 0, . . . , m − 2 另外，弧線函數也有非常廣泛的應用。例如，Runge[20] 發現，使用高階多項式進行插值的時候，在邊界點上會出現劇烈的波動現象，從而導致估計誤差的增大。之後，人們發現，利用弧線函數進行插值處理，可以避免上述現象，從而使得弧線函數在函數近似問題上得到了很好的應用。 1.

(10) “在計算機科學的計算機輔助設計和計算機圖形學中，弧線通常是指分段定義的多項式參數曲線。由於弧線構造簡單，使用方便，配適準確，並能近似曲線配適和交互式曲線設計中複雜的形狀，弧線是這些領域中曲線的常用表示方法1 。” 由於弧線函數在函數近似方面的良好性質，見 de Boor[4]，它在統計學領域的發展也日益深入。而對於如何更準確地估計未知函數則是統計學領域的研究熱點。因此，通過研究比較基於 Splines 函數的不同估計方法則顯得非常有意義。. 1.2 研究目的 [ ] 假設有一個反應變數 Y 和一個在 a, b 之間取值的變量 X，並且進一步假設他們之間存在如下函數關係：. 立. 治政 Y = f (X) + ε 大. (1.2). ‧ 國. 學. 其中，ε 是一個均值為 0，且變異數 E(ε2 ) = σ2 > 0 的隨機變數。通過既有的數據，即 ( ) X, Y 的觀測值，我們可以對 X 和 Y 之間的函數關係 f 進行估計。估計函數的方法有很多，一種是假定待估函數具有某種形式，如估計密度函數時候假定常. ‧. 態分佈，然後通過最大概似或者是動差估計的方法將函數的待估參數估計出來，從而估計函數。另外一種則是將待估計函數透過一連串基底函數表示，如多項式和本文將要討論的弧線函. y. Nat. sit. 數 (spline)。多項式的缺點在於無法很好地配適函數在局部定義域內的形狀，而弧線函數每一個. er. io. 基底函數都只定義局部函數特徵，因此能較好地配適函數的局部特徵，從而得到更精確的 fˆ。. al. v. n. 根據前面一章所給出弧線函數，通過 de Boor [4] 給出的相關算法，我們可以將弧線函數. i n Ch 由以下 m + k 個 B-Spline 的基底函數表達出來： engchi U (

(11) ) (

(12). ) N ·

(13) m, z1 , . . . , zm+1 , . . . , N ·

(14) m, zm+k , . . . , z2m+k ,. 其中，z1 , . . . , z2m+k 是由 a, . . . , a, γ1 , . . . , γk , b, . . . , b 組成的節點序列，而每一個基底函數可以 | {z } | {z } m times. m times. 透過遞迴的方式定義，其定義方式如下： N ( x |m, z1 , . . . , zm+1 ). =. 1. x − z1 z m +1 − x N ( x |m − 1, z1 , . . . , zm ) + N ( x |m − 1, z2 , . . . , zm+1 ) z m − z1 z m +1 − z 2  1, if x ∈ [z1 , z2 ); N ( x |1, z1 , z2 ) = 0, otherwise.. https://zh.wikipedia.org/wiki/弧線函數. 2.

(15)   ( z 2 − x ) m −1 / ( z 2 − z 1 ) m −1 ,. if x ∈ [z1 , z2 );. N ( x |m, z1 , . . . , z1 , z2 ) = | {z } 0, m times   ( x − z 1 ) m −1 / ( z 2 − z 1 ) m −1 , N ( x |m, z1 , z2 , . . . , z2 ) = | {z } 0,. otherwise. if x ∈ [z1 , z2 ); otherwise.. m times. 然後，根據這些基底函數，在 [ a, b] 上的函數我們可以通過 spline basis functions 近似，即得到如下迴歸模型： (

(16) ) (

(17) ) Y = c1 · N X

(18) m, z1 , . . . , zm+1 + · · · + cm+k · N X

(19) m, zm+k , . . . , z2m+k + ε. (1.3). 以上方程有三組未知數需要估計，一是基底函數相對應的迴歸係數，即 c1 , . . . , cm+k , 二是基底函數的 order，三是內部節點的數量與位置，即 γ1 , . . . , γk 。因為迴歸係數的估計相對簡單. 政治大. 因此本文對此不作探討，另外，本文中我們將會把基底函數的 order 固定。因此，本文主要通. 立. 過對比不同的方法在內部節點的數量與位置估計上所導致函數估計精度的差異，進而探討不同. ‧. ‧ 國. 1.3 研究框架. 學. 方法的相對適用環境。. 本文一共包含 5 個章節，其中，第一章為前言，包括研究動機、研究目的與研究框架，勾. sit. y. Nat. 勒整篇文章的大致輪廓。第二章為文獻探討，主要討論與本文相關領域的學術發展情況。第三. io. er. 章為研究方法，解釋本文所要比較的三種方法。第四章為模擬資料分析，分為資料產生與實證結果，分析三種方法各自的表現情況。第五章為結論與建議，根據前面的數據結論進行總結分. al. n 析並給出相對應的建議。. Ch. engchi. 3. i n U. v.

(20) 2 文獻探討使用 spline 近似迴歸函數的時候，節點可以先行固定再對迴歸係數進行估計，例如 Friedman and Silverman[9],Kooperberg et al. [13], Stone et al. [24] 等人在文中將內部節點的位置固定的做法。或是使用 free knot spline，即將節點的位置視為未知參數從而在迴歸分析中將其與迴歸係數一起估計。Burchard[3] 認為使用 free knot spline 能夠提高弧線函數的估計結果。相較於其他一些傳统的弧線函數，如 Eubank[7], Wahba[27]. 或者一些通過增添懲罰項的處理方式，如 O’Sullivan[18], Eilers and Marx[6]. 使用 free knot spline 是一種局部式的調整方式，其内部節點通過移动到函数突變程度較大的地方從而更好地估計原有函數。另外，Smith et al.[23],Denison et al. [5] 等人的貝氏無母數迴歸與使用 free knot spline 的方法較為接近，其. 政治大另一方面，Jupp[10], Lindstrom[15] 立認為將節點視為參數雖然能提高函數的估計精度但卻. 中節點的位置亦為未知參數，其先驗分配值在若干位置不為 0.. ‧ 國. 學. 會帶來兩個主要問題，一是計算最優化的節點的位置會隨著節點數量的增長而越來越困難，二是最優化的節點中往往某些節點重複出現，從而影響估計精度。為了解決這些問題，很多學者提出了不同的方法。. ‧. Lindstrom[15] 提出了一種修正最小平方法的方法，既提升了節點位置的計算速度也保留. sit. y. Nat. 了 free knot spline 的估計精度，從而提高函數的估計效果。Molinari et al.[17] 則將節點的位置限制在某部份區域內從而降低了計算量。Beliakov[1] 將 cutting angle deterministic global 最優. io. n. al. er. 化方法應用到了 spline 的估計問題上，並認為在節點數相對較小的情況下該方法的計算速度較快。. Ch. engchi. i n U. v. 另一部分學者發展了一係列逐步增加或刪除節點的方法。Smith[23] 將變數選取的方法應用到了節點的選取上面，逐步刪除過多的節點。遵循變數選取的思路，Friedman and Silverman[9], Friedman[8] 提出了 TURBO 和 MARS 兩種選取節點的辦法。同樣的還有 Stone et al. [24] 的 POLYMARS，Denisonet al. [5] 的 Bayesian MARS。但 Zhou and Shen[29], Miyata and Shen[16] 認為這些方法只能在有限的範圍內選取內部節點，因此不可能是最優化的節點組合。另外，Zhou and Shen[29] 提出了一種新的選擇最佳節點組合的方法，這種方法對函數突變程度更大的區域賦予更多的節點，與之前逐步增加或者減少的方法相比，這種新方法的誤差更小。 Kang et al.[12] 提出了一個兩步選擇法。首先，通過求解一個 sparse 優化問題得到一個包含較多節點的向量。然後，通過刪除多餘的節點並調整剩下節點的位置從而得到最終的結果。. 4.

(21) 與其他方法相比，作者通過模擬認為他的方法能夠得到較少的節點。但是他的方法計算量較大，且在數據的擾動項較大的情況下表現不佳。與 Kang 的方法採用類似思路的有，Van Loock et al.[25] 提出通過求解 sparse 優化問題從而刪除多餘節點的方法。Yuan et al.[28] 從初始節點集合中選取若干節點作為最終結果。另外，一些別的方法也被提出，如 Lee[14] 通過將最小距離法則 (MDL) 應用在節點數量及位置的估計中。此外，一些學者如 Smith and Kohn[22], Denison et al.[5], Biller[2] 在廣義線性迴歸模型的條件下提出了基於馬爾科夫蒙特卡羅方法的貝氏適應性 spline，從而為 free knot spline 的研究提供了另外一種思路。Kaishev et al.[11] 基於幾何控制多邊形的想法提出了一種新的關於節點的選擇方法，這種新的方法避免了最優化的問題，使得計算效率大大提升。. 政治大. 本文的目標在於透過比較 Lindstrom[15] 與 Kaishev et al.[11] 所提出的兩種不同的方法，從而探討兩種方法彼此的優缺點。. 立. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 5. i n U. v.

(22) 3 研究方法假設對於（1.2）中的 ( X, Y )，我們有 n 筆觀察值，即 ( x1 , y1 ), . . . , ( xn , yn )，並且我們用 B-Spline 基底的線性組合去近似 (1.2) 的迴歸函數 f ，即考慮（1.3）的模型。本章將會討論近似過程中所採用的不同的節點選擇方法。第一節介紹基於最優化思路的節點選擇方法。第二節介紹基於幾何控制多邊形特性的求取節點的方法。第三節為了豐富比較，介紹了固定節點的方法。. 3.1 最佳化方法 3.1.1. 基於 GCV 和修正後 GCV 的估計. 政治大介紹修正後的 GCV 估計節點的方法。當我們用前文所介紹的由基底函數所構成的迴歸方程去立. 本節中我們會先介紹基於 GCV (Generalized Cross Validation) 估計節點的方法，然後再. 近似原函數 f 的時候，若給定內部節點向量 γ，則模型（1.3）簡化成一般線性模型。令 cˆ(γ). ‧ 國. 學. 為給定內部節點向量之後所估計出來的迴歸係數而 R(γ, cˆ(γ)) 為迴歸方程的殘差平方和，則當節點數目給定時，節點向量可以透過最小化目標函數 R(γ, cˆ(γ)) 而解得。. ‧. 當需要決定具體節點數的時候，繼續最小化 R(γ, cˆ(γ)) 這個目標函數容易出現過度配適的. y. Nat. 現象，因此需要對模型的複雜程度進行懲罰。Eubank [7] 提出了一個稱為 GCV 的方法用於選. er. io. sit. 擇節點的位置與個數，其標準等同於最小化. Gss (γ) = R(γ, cˆ(γ))/[trace( I − Ac|γ )]2 = R(γ, cˆ(γ))/(n − nb )2 ,. n. al. Ch. n engchi U. iv. (3.1). 其中，n 為觀察資料的筆數，nb 為基底函數的個數，I 為單位矩陣，Ac|γ 是給定 γ 下迴歸模型 (1.3) 的 influence 矩陣 (投影矩陣)。上述最小化 Gss 得到節點估計的方法就是基於 GCV 的節點估計方法。以下我們再介紹一個修正後的 GCV 的做法。 (3.1) 中的 influence 矩陣是在節點向量 γ 給定之下得到的，因此，當 γ 也需要被估計的時候，Lindstrom [15] 認為需要對 Ac|γ 進行修正。修正後的 Gss (γ) 等同於 G f ull (γ) = R(γ, cˆ(γ))/(n − nb − k)2 .. (3.2). 以下我們簡述一下他的推導過程。(1.3) 的迴歸模型可以寫成如下形式 Y = η (θ ) + ε，θ =. [γT , cT ]T ，則此時的 influence 矩陣 Aθˆ 滿足 η (θˆ) = Aθˆ y, 因此可以寫成如下形式： Aθˆ =. ∂η (θˆ) ∂η (θˆ) ∂θˆ = ∂y ∂θˆ ∂y 6. (3.3).

(23) 其中，y 是 Y 的觀察值向量。θˆ 是 θ 的最小平方估計。由於在非線性的條件下，無法得出具體的矩陣 Aθˆ ，因此需要一個近似解。pregibon[19] 提出了一個解決方案，假設 δ(θ, y) =. ( XθT Xθ )−1 XθT (y − η (θ ))，其中 Xθ = ∂η (θ )/∂θ。給定一個在 θˆ 附近的 θ0 點，θˆ ≈ θ0 + δ(θ0 , y)，並且， ∂θˆ ∂[θ0 + δ(θ0 , y)] ≈ = ( X0T X0 )−1 X0T , (3.4) ∂y ∂y

(24)

(25)

(26) = Xθ0 . 進一步地，如果令 Xˆ = ∂η (θ )/∂θ

(27) ˆ = X ˆ ，則 influence 矩陣可以定義其中 X0 = ∂η θ ∂θ θ θ 0. 成如下： Aθˆ = lim Xˆ ( X0T X0 )−1 X0T = Xˆ ( Xˆ T Xˆ )−1 Xˆ T .. (3.5). θ0 →θˆ. [ ] 利用上述表達式，在模型 (1.3) 中，Xˆ = Xˆ 1 , Xˆ 2 ，其中，Xˆ 1 是 B-Spline 基底函數在解釋變. 政治大 n × k 的矩陣。當 Xˆ 是一個滿秩的矩陣時，立. ˆ Xˆ 2 = (∂ Xˆ 1 /∂γˆ )cˆ(γˆ ) 是一個數的觀察值處取值的 n × (m + k) 矩陣，而基底函數的節點為 γ，. ‧ 國. 學. trace( I − Aθˆ ) = trace( I − Xˆ ( Xˆ T Xˆ )−1 Xˆ T ) = n − nb − k. 因此，(3.1) 的 Gss 可以用 (3.2) 的 G f ull 取代。上述最小化 G f ull 得到節點估計的方法就是基於. ‧. 修正後 GCV 的節點估計方法。. Nat. y. 帶懲罰項的估計. sit. 3.1.2. n. al. er. io. Lindstrom[15] 認為前述 G f ull 的估計方法並沒有處理內部節點可能過近的問題，因此提出. v. 了一個帶有懲罰項且可提高計算效率的目標函數，藉由最小化此目標函數得到節點估計，其目標函數如下：. Ch. engchi. i n U. R(γ, cˆ(γ)) · J (γ, p),. (3.6). 其中，p 是用以控制懲罰項 J (γ, p) 大小的參數， J (γ, p) = α( p, γ0 (k)) · log( P(γ)) + 1, α( p, γ0 (k )) = ( p − 1)/[log( P[γ0 (k )])], k +1. log( P(γ)) = − ∑ log((k + 1) · hi ), i =1. h i = ( γi − γi − 1 ) / ( b − a ). i = 1, ..., k + 1,. γ0 = a, γk+1 = b, (γ1 , . . . , γk ) = γ.. 7. (3.7).

(28) γ0 (k ) 則是一組固定的內部節點向量，長度為 k，它所對應的 h1 , . . . , hk+1 如下：h1 = d/(k + 1)，hi − hi−1 = [1 − d/(k + 1)]/k, i = 2, . . . , k + 1. 本文中 d = 0.04, a = 0, b = 1. 如果節點非常接近，則公式 3.7中的 hi 會變得很小從而 log( P(γ)) 變很大，最終導致懲罰項 J (γ, p) 變很大。另外，Lindstrom 認為由於懲罰項的引入，上一節中所推導出的 Aθˆ 無法體現自由度的變 [ pen]. pen. 化，因此提出了一個新的選擇標準 GLS (θ ) = R(γ, c)/[trace( I − Aθˆ )]。以下將定義 R(γ, c) pen. 以及 Aθˆ . 對於 γ = (γ1 , . . . , γk ), c = (c1 , . . . , cm+k ), x ∈ [ a, b], 令 f ( x |c, γ) = ∑ik=+1m ci Bi ( x ), 其中 B1 , . . . , Bm+k 為在 [ a, b] 上且內部節點向量為 γ 的 B-Spline 基底函數。則 R(γ, c) = ∑in=1 (yi − pen fˆ( xi |c, γ))2 , 其中 ( x1 , y1 ), . . . , ( xn , yn ) 是 ( X, Y ) 的觀察值。Aθˆ 具體推導過程如下：. 治政 Gradient) 和二次偏微 (即 Hessian). 因此，公式 (3.4) 會變成：大立 ∂θˆ ∂(θ + δ(θ , y)) ∂Ω. 令 δ(θ, y) = −Ω−1 w，其中 w 和 Ω 分別是新的目標函數 (3.6) 關於 θ 的一次偏微 (即. ‧ 國. 0. ∂y. =−. −1. ∂y. |θ0 w0 − Ω0−1. ∂w |θ , ∂y 0. 學. ∂y. 0. ≈. pen. 其中，Ω0 與 w0 分別是 Ω 與 w 在點 θ0 上的取值。然後可以得出新準則中的 Aθˆ ： ( ) ∂Ω−1 ∂w ˆ −1 · ∂w |θˆ , = lim Xˆ − |θ0 w0 − Ω0−1 |θ0 = − Xˆ · Ω ∂y ∂y ∂y θ0 →θˆ. ‧. pen. Aθˆ. Nat. sit. y. ˆ 是 Ω 在點 θˆ 的取值。令 r = y − η (θˆ),J = J (γ, ˆ p ), 則其中，Ω. n. al. er. io. ∂w | ˆ = 2(− J · Xˆ T + ∂J (γ, p)/∂θ |γ=γˆ · r T ), ∂y θ. i C ˆ = 2 · J · Xˆ Xˆ + 0 0 , n Ω h e n g c h0 i S U [. ]. v. T. S = −2( Xˆ 2T · r · ∂J/∂γT |γ=γˆ + ∂J/∂γT |γ=γˆ · r T · Xˆ 2 ) + r T r · ∂2 J/∂γ2 |γ=γˆ . 由於參數 p 可以控制懲罰的大小，所以 p 的選擇會影響到目標函數的選擇。當 p ≤ 1.1 的時候，Lindstrom 認為應該使用無懲罰項的目標函數並且以 G f ull 作為選擇標準。當 p > 1.1 的 [ pen]. 時候，應該使用帶有懲罰項的目標函數並且以 GLS (θ ) 作為選擇標準，即： Gcomb (θ, p) =.   G f ull (θ ),. if p ≤ 1.1;.  G [ pen] (θ ), if p > 1.1. LS. 通過以上目標函數以及判斷標準就可以選擇出符合條件的節點數量及其位置。. 8.

(29) 3.2. GeD 方法 kaishev et al. [11] 提出 GeD(Geometrically designed spline regression method) 方法，根. 據資料自動決定節點的位置與階數的大小。以下我們假設 X 的觀察值已經經過排序，即 x1 ≤ x2 ≤ · · · ≤ xn . 其次，由於本文統一在 Cubic-Spline 的架構下進行比較，因此，將不進行階數的調整，並設定階數等於 4。其具體算法如下： Step 1. 考慮內部節點數為 0 的情形。此時，B-Spline 的基底函數有 4 個，先由現有數據配適出迴歸函數的估計 fˆ，接著計算出殘差 ri = yi − fˆ( xi ), i = 1, . . . , n, 和殘差平方和 ∑in=1 ri2 . Step 2. 根據每一個殘差的正負符號將其與鄰近殘差進行分群。令 u 為群數，d j 為第 j 群中的殘差個數，j = 1, . . . , u, 即：. 政治大立) =sign(r ) = · · · = sign(r ). sign(r1 ) = · · · = sign(rd1 ). ̸= sign(rd1 +1. d1 +2. d1 + d2. ‧ 國. 學. .. .. ̸= sign(rd1 +d2 +···+du−1 +1 ) =sign(rd1 +d2 +···+du−1 +2 ) = · · · = sign(rd1 +d2 +···+du−1 +du ). ‧. 其中，∑uj=1 d j = n.. Nat. xd( j) − xd( j−1)+1 , d( j) = d1 + d2 + · · · + d j ,. y. Step 3. 分別計算這 u 個殘差群的群內均值 m j 以及殘差對應的 X 觀察值的範圍，即 η j =. sit. j = 1, . . . , u.. al. er. io. Step 4. 通過 mmax = max1≤ j≤u (|m j |) 和 ηmax = max1≤ j≤u (η j ) 將各個群內均值的絕對值和. n. 群範圍標準化得到 m′j = |m j |/mmax , 0 < m′j ≤ 1 和 η j′ = η j /ηmax , 0 < η j′ ≤ 1. Ch. Step 5. 計算各群的權重 w j = β · m′j + (1 − β) · η j′ ,. engchi. i n U. v. j = 1, . . . , u. 這個權重是用來衡量各個. 群內的迴歸函數估計值和真實 f 值之間偏離情況的指標，其中，β 是模型參數並事先給定，在本文中，它被設定為 0.5。 Step 6. 對權重 w1 , w2 , w3 , . . . , wu 進行排序，假設最大的權重為第 jρ 群的權重，然後檢查 [ ] 第 jρ 群殘差所對應的 X 的觀察值的範圍是否已經存在節點，即是否存在 γi ∈ xd( jρ −1) , xd( jρ −1)+d jρ ， [ ] 如果沒有，就在 xd( jρ −1) , xd( jρ −1)+d jρ 中插入一個新的節點 γ∗ =. ( d( jρ −1)+d jρ. ∑. i =d( jρ −1). ri xi. )/( d( jρ −1)+d jρ ) ∑ ri . i =d( jρ −1). [ ] 如果在區間 xd( jρ −1) , xd( jρ −1)+d jρ 已經存在節點，則找到第二大的權重，然後重複上述步驟。 Step 7. 通過新的節點向量配適出一個新的估計函數，並且計算 RSS，即殘差平方和。. 9.

(30) Step 8. 令 RSS(k ) 為 k 個節點下配適出來的迴歸模型的殘差平方和，計算如下指標： π=. RSS(目前節點數) RSS(目前節點數 − q). 其中，q 是事先給定的大於或等於 1 的整數。當 π 小於給定的判斷們檻值 πexit 時，則回到 Step 2，反之停止插入新的節點。在本文中 q = 2，πexit = 0.9。另外，如果 k ≤ q，那麼整個運算回到 Step 2. 再循環一次。通過上述算法，最後我們可以求得具體需要的節點數量與相對應的位置。. 3.3 固定等距節點. 政治大. 本文也考慮固定節點的基礎方法作為與之前兩種使用 free knot 的方法進行對比。 [ ] 假設給定 k − 1 個節點作為配適弧線迴歸函數的內部節點，且在區間 a, b 上進行插入，. 立. 那麼每個具體節點的位置將會如下：. ‧ 國. 學. b−a k b−a γ2 = a + 2 · k .. . γ1 = a +. ‧ y. b−a k. sit. Nat. γk = a + ( k − 1 ) ·. al. er. io. 然後，通過 BIC(Bayesian Information Criterion) 對不同個數的節點組合進行挑選，對應 BIC. v. n. 最小值的那一組節點將會成為最後的內部節點向量。其中，BIC 的具體公式如下：. i n C h=elnn(ng BIC )k − h c 2lni(Lˆ )U 0. Lˆ 是根據觀測資料所能求出的概似函數的最大值，k0 為模型中被估計的參數個數，即等於 k + 3.. 10.

(31) 4 模擬資料分析 4.1 資料產生與模擬流程首先，所有資料均由模型 (1.3) 生成，具體參數、步驟如下： [ ] Step 1.X 的觀察值通過在 0, 1 中等距劃分方式產生。共有兩組，一組較為稀鬆，為 40 筆，一組較為緊密，為 400 筆。用以驗證資料的多寡程度對於估計方法的影響。 ( ) Step 2. 利用 Uni f orm 0, 1 隨機產生 KG 個點作為真實資料的內部節點。在本文中，這樣的節點數一共有 3 組，分別是 3 個、5 個和 7 個。用以驗證函數突變程度對於估計方法的影響。 ( ) Step 3. 利用 N 0, 1 隨機產生模型 (1.3) 中的迴歸係數，並且對這些係數進行標準化處理使與其相對應的迴歸函數值介於 0, 1 之間。在本文中，這樣的迴歸係數一共有 3 組，分別對應. 政治大. 節點數 3、5、7 的組別。 ( Step 4. 利用 N 0, σ2 ) 隨機生成擾動項 ε，在本文中，σ 會取 0.1 和 0.01 用以驗證誤差大. 立. ‧ 國. 學. 小對於估計方法的影響。. Step 5. 根據第一步中產生資料和第二步中的內部節點，我們可以計算出相對應的 B-Spline. 的模擬資料。. ‧. 基底函數，再乘以第三步中的迴歸係數和加上第四步中的擾動項就可以得到本文中需要使用到. Nat. sit. y. 其次，本文將會以上述模擬資料驗證前文論述的三種方法. io. al. n. 下：. er. 1. 對於帶懲罰項的估計方法，根據 p = 1.1, 2, 11, 101 來分別進行模擬運算，具體過程如. i n U. v. Step a. 針對內部節點數 k = 2, . . . , 10，分別記錄能夠最小化 R(θ ) J (θ, p) 的節點向量。在. Ch. engchi. 此過程中，由於多元函數求極值存在一些不穩定性，因此每一次求最值的運算都會重複兩次，然後取較小的那次作為最後的結果。 Step b. 將不同 k 下最優化的節點向量代入相對應得評判標準中，記錄下結果最小的那一組節點向量。 Step c. 重複前面兩步 10 次。 2. 對於 GeD 方法，則根據 πexit = 0.9 重複 10 次模擬計算並記錄模擬結果。 3. 對於固定等距節點方法，由於節點位置固定，直接代入資料，重複 10 次並記錄模擬結果。. 11.

(32) 4.2 實證結果 4.2.1. 節點選擇比較. Case 1. 資料較少且誤差較大的情況我們首先對各個方法估計真實節點數的準確程度進行比較分析。圖 (4.1),(4.2),(4.3) 是在資料數較少且誤差較大的情況下各方法對真實節點數做出的估計：. 立. 政治大. ‧. ‧ 國. 學. n. al. er. io. sit. y. Nat. 圖 4.1: σ = 0.1,n = 40,KG = 3 下節點估計數. Ch. engchi. i n U. v. 圖 4.2: σ = 0.1,n = 40,KG = 5 節點估計數. 12.

(33) 政治大. 圖 4.3: σ = 0.1,n = 40,KG = 7 下節點估計數. 立. ‧ 國. 學. 從上述圖片中可以發現，在資料筆數較少的情況下，無論真實節點數如何變化，每一種方法的所估計的節點數趨勢較為一致，即 G f ull 傾向於低估節點數，而其他方法則傾向於高估。但在真實節點數較多的，即 K g = 7 的情況下，使用固定節點位置和 BIC 準則的方法估計的節點. ‧. 數較為接近真實節點數。這可能是因為缺乏足夠的資料去進行函數的最優化運算從而導致了其. io. y. sit. 資料確定節點位置。. Nat. 他方法的傾向於通過高估節點數來降低迴歸方程的最小平法和，而固定位置的方法不需要依賴. n. al. er. Case 2. 資料較多但誤差較大的情況. i n U. v. 維持誤差不變，當我們調整資料筆數，從 40 筆增加到 400 筆的時候，從圖 (4.4),(4.5),(4.6). Ch. engchi. 中我們可以看到節點估計的情況出現了一些改變：. 13.

(34) 政治大. 圖 4.4: σ = 0.1,n = 400,KG = 3 下節點估計數. 立. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. 圖 4.5: σ = 0.1,n = 400,KG = 5 下節點估計數. 14.

(35) 政治大. 圖 4.6: σ = 0.1,n = 400,KG = 7 下節點估計數. 立. ‧ 國. pen. 是多少，其節點數量估計表現較為準確。而 Gls. 學. 從圖中可以看到，在資料數量增加之後，Geometrically Designed 的方法無論真實節點數在真實節點數較少的情況下容易高估節點數. 量，這點與 G f ull 相反。此外，在真實節點數較高的情形下，即 K g = 7，除固定位置的方法外，. ‧. 其餘的方法的表現均與正確答案較為接近。. sit. y. Nat. Case 3. 資料較多且誤差較小的情況. io. 反應了估計情況的改變：. n. al. er. 我們進一步調整資料，維持資料筆數不變，將資料的誤差縮小到 0.01，圖 (4.7),(4.8),(4.9). Ch. engchi. i n U. v. 圖 4.7: σ = 0.01,n = 400,KG = 3 下節點估計數. 15.

(36) 政治大. 圖 4.8: σ = 0.01,n = 400,KG = 5 下節點估計數. 立. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. 圖 4.9: σ = 0.01,n = 400,KG = 7 下節點估計數. 從圖中可以看到，誤差的減小大幅降低了 Geometrically Designed 方法和固定節點位置方法對於節點數量的估計準確度，無論真實節點數是多少，這兩種方法都會過高地估計節點數 pen. pen. 量。另外，除了在 K g = 5 的情形下 Gls 相對高估了節點數量外，Gls 和 G f ull 都能非常接近地估計正確節點的數量。. 16.

(37) 4.2.2. 估計精度比較. 除了節點數量的估計，對於函數估計來說，更為重要的應該是估計出來的函數與原函數相差多少，下面，我們將針對不同情況對函數的估計精度進行比較。其中，Mean of Loss 是指 Mean Squared Error，Time User 是指程式運行所花費時間。 Case 1. 資料較少且誤差較大的情況針對資料數量較少的情況，表 (4.1),(4.2),(4.3) 呈現了各種方法的估計情況：表 4.1: k G = 3, σ = 0.1, n = 40 σ = 0.1, k G = 3, fixed coefficients and knots Geometrically Designed. pen. pen. Gls. Gls. p=2 政治 9 10大 10 8 8 9 10 3 8 5 αexit = 0.9. 立99565. 9 10 10 8 8. 99997. 9 10 10 9 7. (第 4,7,8 次異常). (第 3 次異常). 學. ‧ 國. Number of Knots. p=3. 0.3092914. 0.3083854. Mean of Loss. 0.005589562. 0.08498319. 0.009034612. Time User. 0.577. 937.068. 908.714. BIC with. pen Gls. Equally Spacing. p=101. Time User. n. Mean of Loss. al. y. sit. (第 1,3 次異常). er. io Mean of R(θ ). 10 10 9 9 7. 55433. 0.3225177 iv n C 0.003550839 U h e n g c h i 0.01575261 0.392 992.500 0.3170897. 17. G f ul p=1.1. 9 10 10 8 8. 63357. Number of Knots. ‧. 0.3002208. Nat. Mean of R(θ ). 22224 22222 0.299624 0.00262018 132.099.

(38) 表 4.2: k G = 5, σ = 0.1, n = 40 σ = 0.1, k G = 5, fixed coefficients and knots pen. Geometrically Designed. Gls. p=2. p=3. 8 8 10 10 7. 8 10 10 10 9. 10 10 9 10 10. 10 10 10 10 9. (第 6,8,10 次異常). (第 6 次異常). αexit = 0.9 8 18 14 12 11. Number of Knots. 8 7 4 12 11. pen. Gls. Mean of R(θ ). 0.2490128. 0.2335129. 0.2313402. Mean of Loss. 0.006776192. 0.1060095. 0.007291692. Time User. 0.923. 923.925. 905.590. BIC with. pen Gls. G f ul. Equally Spacing. p=101. p=1.1. 政治 9 10大10 10 9 10 18 12 12 7 10 10 10 10 9 立5 5 7 7 10 (第 6 次異常). Number of Knots. 33463 33334. 0.2304437. 0.2522541. Mean of Loss. 0.01525978. 0.004943531. 0.003202004. Time User. 0.383. 892.237. 119.333. ‧. ‧ 國. 0.2425759. 學. Mean of R(θ ). n. er. io. sit. y. Nat. al. Ch. engchi. 18. i n U. v.

(39) 表 4.3: k G = 7, σ = 0.1, n = 40 σ = 0.1, k G = 7, fixed coefficients and knots pen. Geometrically Designed αexit = 0.9. Gls. p=2. p=3. 10 5 10 10 10. 9 13 10 9 11. Number of Knots. pen. Gls. 10 5 10 10 10. 10 9 6 10 10. 5 10 5 5 6. (第 9 10 次異常). 7 9 10 8 9. Mean of R(θ ). 0.2906263. 0.3319193. 0.327041. Mean of Loss. 0.0199668. 0.01080632. 0.01196568. Time User. 0.722. 984.259. 959.288. BIC with. pen Gls. G f ul. Equally Spacing. p=101. p=1.1. 10 10 10 政治 10 5大 9 4 10 13 3 7 9 10 8 10 44443. 52226. Mean of R(θ ). 0.3646743. 0.3276086. 0.3677975. Mean of Loss. 0.01028184. 0.01013793. 0.00628123. Time User. 0.424. 885.792. 95.357. Number of Knots. 22222. (第 10 次異常). 學. ‧. ‧ 國. 立. pen. 結果可以看到，在資料數量較少的情形下，Gls. io. sit. y. Nat. 表中標示所謂異常狀況是指當次估計的 MSE 非常大，即與其他結果差異較大的情況。從的估計結果非常不穩定，且在真實節點數較. al. n. pen. 好。對於 Gls. er. 少的情形下估計精度不如 G f ull 和 Geometrically Designed 的方法，其中又以 G f ull 的表現最. i n U. v. 出現的異常情況，我們以 σ = 0.1, KG = 3, p = 2, n = 40 的第四次異常情況和. Ch. engchi. σ = 0.1, KG = 3, p = 2, n = 400 的正常情況作對比2 ，如圖 (4.10),(4.11)： 2. 兩次對比的模擬資料與運算過程均通過 R 語言內的 set.seed 進行固定. 19.

(40) 政治大. 圖 4.10: 異常情況下估計表現. 立. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. 圖 4.11: 正常情況下估計表現. 可以看到，隨著數據量的增長，函數估計的精確度有了大幅提高，相較於由於數據太少導致估計偏差的情況，n = 400 下的估計出的函數已經基本與原函數重合。此外，對三種方法所需時間，由於 Geometrically Designed 的方法不需要進行非線性函數 pen. 的優化問題，因此運行時間遠遠小於其他方法。而 Gls. 並沒有表現出其應有的在計算速度方. 面的提升。 Case 2. 資料較多但誤差較大的情況針對資料過少的問題，我們通過將資料數增加到 400 筆，進而比較各方法估計效果的改 20.

(41) 變，結果如表 (4.4),(4.5),(4.6)：表 4.4: k G = 3, σ = 0.1, n = 400 σ = 0.1, k G = 3, fixed coefficients and knots Geometrically Designed. pen. pen. Gls. Gls. αexit = 0.9. p=2. p=3. 65555. 9 8 2 2 10. 9 10 10 8 10. 53555. 10 10 2 9 8. 2 9 3 10 10. Mean of R(θ ). 4.15106. 4.031897. 4.099653. Mean of Loss. 0.0008639218. 0.000711288. 0.0009825077. Time User. 1.678. 8164.576. 7947.617. BIC with. pen Gls. G f ul. Number of Knots. p = 101 治政 65658 9 10 10 10 10 大 66666 9 9 3 10 10. 32225. 3.890481. 4.129813. 3.876426. Mean of Loss. 0.0003044727. 0.001100262. 0.0002577466. Time User. 0.533. 8392.959. 261.559. Equally Spacing Number of Knots. 22223. 學. al. αexit = 0.9. pen. Gls. sit. io. Geometrically Designed. er. Nat. 表 4.5: k G = 5, σ = 0.1, n = 400 σ = 0.1, k G = 5, fixed coefficients and knots. y. ‧. ‧ 國. Mean of R(θ ). 立. p = 1.1. p=2. pen. Gls. p=3. Mean of Loss. 0.0013043768. 0.0003933684. 0.0003817892. Time User. 1.678. 8729.252. 8312.817. BIC with. pen Gls. G f ul. Equally Spacing. p = 101. p = 1.1. 77777. 9 10 10 10 7. 43338. 7 7 7 10 7. 10 9 9 9 10. 33543. Mean of R(θ ). 4.11438. 3.902046. 3.922526. Mean of Loss. 0.0007935589. 0.000392659. 0.0003685712. Time User. 0.590. 8568.156. 200.480. n Mean of R(θ ). 67676 9 10 9 10 i v10 n C4h6 5 4 8 10U 10 9 9 10 i e h n c g 4.37938 3.903166. Number of Knots. Number of Knots. 21. 9 10 10 10 10 10 10 9 9 10 3.908453.

(42) 表 4.6: k G = 7, σ = 0.1, n = 400 σ = 0.1, k G = 7, fixed coefficients and knots Geometrically Designed. pen. pen. Gls. Gls. αexit = 0.9. p=2. p=3. 56588. 8 10 4 6 10. 8 8 4 6 10. 75879. 10 6 10 9 10. 10 9 10 9 10. Mean of R(θ ). 4.332343. 4.929553. 4.943285. Mean of Loss. 0.001460695. 0.002827422. 0.002861437. Time User. 2.248. 7225.801. 7445.85. BIC with. pen Gls. G f ul. Equally Spacing. p = 101. p = 1.1. 10 14 10 10 10. 77359. 8 6 4 5 10. Number of Knots. Number of Knots Mean of R(θ ). 立. Mean of Loss. 治 98989 政大 4.077527 4.943285. 13 13 10 13 10. ‧ 國. 3.99468. 0.001030896. 0.002861437. 0.0008321905. 0.652. 7392.553. 209.029. 學. Time User. 69998. pen. 方法不適宜用在資料數較. ‧. 通過增加資料數，估計過程中的異常情況沒有發生，表明 Gls. pen. 少的情況中。此外，Geometrically Designed 和固定節點的方法在這種資料下表現和 Gls. Nat. y. 差. sit. 不多，特別是在 KG = 5 的情況下，後者的精度在 0.0003 左右，遠遠超過 Geometrically De-. al. er. io. signed 的 0.0013 和固定節點的 0.0007。但是，G f ull 在估計精度上仍然是最高的。. n. Case 3. 資料較多且誤差較小的情況. Ch. i n U. v. 由於各種方法本身的估計誤差的絕對值就比較小，因此，如果數據的誤差太大，容易對資. engchi. 料分析產生影響，因此接下來的分析中我們將維持數據量不變，但把數據的誤差縮小至 0.01。結果如表 (4.7),(4.8),(4.9)：. 22.

(43) 表 4.7: k G = 3, σ = 0.01, n = 400 σ = 0.01, k G = 3, fixed coefficients and knots Geometrically Designed. pen. pen. Gls. Gls. αexit = 0.9. p=2. p=3. 11 10 10 10 10. 32222. 33333. 9 10 10 11 10. 22222. 33333. Mean of R(θ ). 0.03864257. 0.04330392. 0.04645359. Mean of Loss. 3.838606e-06. 1.315635e-05. 2.141573e-05. Time User. 2.683. 7392.553. 7594.829. BIC with. pen Gls. G f ul. Equally Spacing. p = 101. p = 1.1. 99999. 23332. 32434. Number of Knots. 治 23332 政大 0.03887681 0.07808387. Number of Knots Mean of R(θ ). 99999. 立 4.117914e-06. Mean of Loss Time User. 0.0001005053. 1.022818e-05. 7146.424. 217.186. ‧. ‧ 國. 0.04096293. 學. 0.538. 93634. 表 4.8: k G = 5, σ = 0.01, n = 400 σ = 0.01, k G = 5, fixed coefficients and knots. n Mean of R(θ ) Mean of Loss Time User. y. pen. Gls. αexit = 0.9. p=2. 9 11 10 10 10. 7 10 9 10 10. 9 10 9 10 10. iv 0.04020865 n C0.04234659 U h e n g c h i 7.37691e-06 1.499823e-05. 8 9 10 10 10. 9 10 10 10 12. sit. io. al. pen. Gls. er. Nat Number of Knots. Geometrically Designed. 8 10 10 10 10. p=3. 0.04105761 1.223566e-05. 2.587. 7190.879. 7186.274. BIC with. pen Gls. G f ul. Equally Spacing. p = 101. p = 1.1. 20 15 15 15 15. 9 10 9 10 10. 54447. 15 15 15 15 15. 8 9 8 10 10. 44545. Mean of R(θ ). 0.04037496. 0.04097288. 0.03913916. Mean of Loss. 9.314688e-06. 1.225837e-05. 3.514883e-06. Time User. 0.617. 7115.406. 203.934. Number of Knots. 23.

(44) 表 4.9: k G = 7, σ = 0.01, n = 400 σ = 0.01, k G = 7, fixed coefficients and knots pen. Geometrically Designed. pen. Gls. Gls. αexit = 0.9. p=2. p=3. 14 15 14 14 15. 8 10 6 9 10. 9 7 6 9 10. 14 14 14 15 15. 10 6 8 10 5. 10 9 7 10 10. Mean of R(θ ). 0.04096927. 0.9989168. 0.9982236. Mean of Loss. 2.115072e-05. 0.002625663. 0.002629872. Time User. 2.587. 7190.879. 7186.274. BIC with. pen Gls. G f ul. Equally Spacing. p = 101. p = 1.1. 21 21 21 21 21. 8 10 6 9 10. 10 9 9 10 10. Number of Knots. Number of Knots Mean of R(θ ). 治 10 9 7 10 8 政大 0.06417372 1.006747. 21 21 21 21 21. 立. Mean of Loss. 8 9 7 10 8 0.04983757. 0.002643602. 6.467929e-05. 0.612. 7191.896. 229.877. ‧ 國. Time User. 學. 9.948728e-05. pen. ‧. 但是，Gls 在 KG = 降低誤差之後，沒有哪種方法的表現能夠在所有情形下超過其他方法。 7 的情形下估計精度較差，與其他方法的差距較大。我們更進一步研究各種方法在沒有誤差的. y. Nat. io. sit. 情形下的估計精度，結果如表 (4.10)：. n. al. Mean of Loss pen Gls , p = 2 pen Gls , p = 11 pen Gls , p = 101. er. 表 4.10: 沒有誤差的情形下各方法表現情況，n=400 K =3 K = 5i v n C h1.096404e-05 5.9648e-05 U i e n g c h0.0003085712 1.030997e-05 G. G. KG = 7 0.002621016 0.002643483. 1.168991e-05. 0.000899794. 0.002629312. G f ull , p = 1.1. 7.387328e-10. 1.447358e-15. 2.173786e-05. Geometrically Designed. 4.429891e-17. 2.064227e-10. 7.002945e-09. 可以看出，無論真實內部節點有多少，Geometrically Designed 的估計精度都比較高。但 pen. 如果資料誤差較大，可能掩蓋掉這部分精度從而導致估計誤差的增加。此外，Gls 的方法在真實內部節點較多的情況下誤差較高。. 24.

(45) 5 結論與建議 5.1 研究結論對於 Lindstrom 提出的 Penalized Estimation 方法，當選擇模型參數 p > 1.1 的時候，即通過最小化 R(γ, cˆ(γ)) · J (γ, p) 求解節點位置，在數據資料較少且誤差較大的情況下，此方法的估計效果不佳，且經常會出現偏差非常大的情況。在數據資料較多但是數據誤差較大的情況下，異常情況不再出現，說明此方法可以在數據較多的情況使用，此外，對於真實內部節點較少的情況，這種方法的表現相較 Geometrically Designed 的方法要來得好，但是，根據模擬結果可以看到，隨著真實內部節點數的增加，此方法的估計誤差隨之增加，因此，可以認為，此. 政治大. 方法不適用於內部節點數較多的資料。在數據資料較多且數據誤差較小的情況下，如果真實內部節點較少，那麼此方法與 Geometrically Designed 在估計準確度上並沒有太大的差異。而對. 立. 於真實內部節點較多的情況，此方法估計誤差仍舊較大。另外，通過模擬研究，隨著 p 的增. ‧ 國. 學. 大，節點個數並沒有隨著減少，估計誤差也沒有呈現出規律性變化。因此可以認為 p 的選擇並沒有起到應有的調節功能。還有一個不足之處就是，這個方法的運算時間非常巨大，特別是. 巨大。. ‧. 10 個最優化節點解的計算特別需要花費時間，另外與之相對應得選擇標準的計算量也是比較. sit. y. Nat. 當選擇模型參數 p ≤ 1.1 的時候，即通過最小化 R(γ, cˆ(γ)) 求解節點位置，最後的結果表. io. er. 明，雖然在內部節點數的估計上，此方法經常會出現低估的現象，但在估計準確度上基本上都要比其他方法要好。此外，這個方法所需要的計算時間並不算太長，比起 p > 1.1 的情況有巨. al. n. 大優勢。. Ch. engchi. i n U. v. 對於 Kang et al. 提出的 Geometrically Designed 方法，從估計誤差上看，相較於 G f ull 、 pen. 和 Equally Spacing，基本排在中間的位置，並沒有太突出的表現。從內部節點數的估計. Gls. 上看，在誤差較小的情況下，Geometrically Designed 會出現高估很多的情況。但是這個方法的優點在於運算時間的大幅降低，由於不需要解非線性函數的最佳化問題，因此計算量大大降低。綜合來看，Penalized Estimation 並沒有表現出文章裡表述的計算效率，相反，運算時間遠遠超過沒有懲罰項的估計方法。另外，從估計誤差上看，帶有懲罰項的估計方法也經常高於沒有懲罰項的估計方法。由於沒有懲罰項的估計方法從最小化目標函數的構建 R(γ, cˆ(γ))，到最後選擇標準 G f ull = R(θ )/(n − nb − k )2 的運用都沒有使用到文中所提及的方法，所以可以認為，文中的方法並不適用於本文中所使用的數據資料。另一方面，Geometrically Designed 確實. 25.

(46) 體現出了其在運算時間上的巨大優勢，因此，如果對於估計誤差要求不高，推薦使用此方法。. 5.2 研究建議 [ ] 本文採用的數據是在區間 0, 1 裡選取固定位置的數據集生成的，因此有可能 Penalized Estimation 的方法在此類數據集中並不適用。因此，可以考慮通過不同類型的數據集進行更細緻化的分析。本文在計算最小化目標函數時使用的是 R Language 內置的 optim 函數，可能由於此函數本身存在一些求解的問題，所以可以通過其他專們的 package 嘗試改善估計結果。本文在模擬運算時採用 R Language 內置的 for 迴圈，而 for 迴圈本身就有速度慢的缺點，. 政治大. 因此可以考慮通過其他計算機語言來進行模擬運算。. 立. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 26. i n U. v.

(47) 參考文獻 [1] Gleb Beliakov. Cutting angle method–a tool for constrained global optimization. Optimization Methods and Software, 19(2):137–151, 2004. [2] Clemens Biller. Adaptive bayesian regression splines in semiparametric generalized linear models. Journal of Computational and Graphical Statistics, 9(1):122–140, 2000. [3] Hermann G Burchard. Splines (with optimal knots) are better. Applicable Analysis, 3(4):309–319, 1974. [4] Carl De Boor. A practical guide to splines; rev. ed. Applied mathematical sciences. Springer, Berlin, 2001.. 立. 政治大. [5] DGT Denison, BK Mallick, and AFM Smith. Automatic bayesian curve fitting. Journal. ‧ 國. 學. of the Royal Statistical Society: Series B (Statistical Methodology), 60(2):333–350, 1998. [6] Paul HC Eilers and Brian D Marx. Flexible smoothing with b-splines and penalties.. ‧. Statistical science, pages 89–102, 1996.. sit. y. Nat. [7] Randall L Eubank. Spline smoothing and nonparametric regression. Marcel Dekker, 1988.. io. pages 1–67, 1991.. er. [8] Jerome H Friedman. Multivariate adaptive regression splines. The annals of statistics,. al. n. iv n C [9] Jerome H Friedman and Bernard W Silverman. Flexible h e n g c h i U parsimonious smoothing and additive modeling. Technometrics, 31(1):3–21, 1989. [10] David LB Jupp. Approximation to data by splines with free knots. SIAM Journal on Numerical Analysis, 15(2):328–343, 1978. [11] Vladimir K Kaishev, Dimitrina S Dimitrova, Steven Haberman, and Richard J Verrall. Geometrically designed, variable knot regression splines. Computational Statistics, 31(3):1079–1105, 2016. [12] Hongmei Kang, Falai Chen, Yusheng Li, Jiansong Deng, and Zhouwang Yang. Knot calculation for spline fitting via sparse optimization. Computer-Aided Design, 58:179–188, 2015. 27.

(48) [13] Charles Kooperberg, Charles J Stone, and Young K Truong. Hazard regression. Journal of the American Statistical Association, 90(429):78–94, 1995. [14] Thomas CM Lee. Regression spline smoothing using the minimum description length principle. Statistics & probability letters, 48(1):71–82, 2000. [15] Mary J Lindstrom. Penalized estimation of free-knot splines. Journal of Computational and Graphical Statistics, 8(2):333–352, 1999. [16] Satoshi Miyata and Xiaotong Shen. Free-knot splines and adaptive knot selection. Journal of the Japan Statistical Society, 35(2):303–324, 2005.. 政治大 regression splines. Computational statistics & data analysis, 45(2):159–178, 2004. 立. [17] Nicolas Molinari, Jean-François Durand, and Robert Sabatier. Bounded optimal knots for. ‧ 國. 學. [18] Finbarr O’sullivan, Brian S Yandell, and William J Raynor Jr. Automatic smoothing of regression functions in generalized linear models. Journal of the American Statistical. ‧. Association, 81(393):96–103, 1986.. [19] Daryl Pregibon. Logistic regression diagnostics. The Annals of Statistics, pages 705–724,. io. sit. y. Nat. 1981.. n. al. er. [20] Carl Runge. Über empirische funktionen und die interpolation zwischen äquidistanten. i n U. v. ordinaten. Zeitschrift für Mathematik und Physik, 46(224-243):20, 1901.. Ch. engchi. [21] Larry Schumaker. Spline functions: basic theory. Cambridge University Press, 2007. [22] Michael Smith and Robert Kohn. Nonparametric regression using bayesian variable selection. Journal of Econometrics, 75(2):317–343, 1996. [23] Michael Smith, Chi-Ming Wong, and Robert Kohn. Additive nonparametric regression with autocorrelated errors. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(2):311–331, 1998. [24] Charles J Stone, Mark H Hansen, Charles Kooperberg, Young K Truong, et al. Polynomial splines and their tensor products in extended linear modeling: 1994 wald memorial lecture. The Annals of Statistics, 25(4):1371–1470, 1997. 28.

(49) [25] Wannes Van Loock, Goele Pipeleers, Joris De Schutter, and Jan Swevers. A convex optimization approach to curve fitting with b-splines. IFAC Proceedings Volumes, 44(1):2290– 2295, 2011. [26] Grace Wahba. A Survey of Some Smoothing Problems and the Method of the Generalized Cross-validation for Solving Them. University of Wisconsin, Department of Statistics, 1976. [27] Grace Wahba. Spline models for observational data. SIAM, 1990. [28] Yuan Yuan, Nan Chen, and Shiyu Zhou. Adaptive b-spline knot selection using multi-. 政治大 [29] Shanggang Zhou and Xiaotong Shen. Spatially adaptive regression splines and accurate 立 resolution basis set. IIE Transactions, 45(12):1263–1277, 2013.. knot selection schemes. Journal of the American Statistical Association, 96(453):247–259,. ‧. ‧ 國. 學. io. sit. y. Nat. n. al. er. 2001.. Ch. engchi. 29. i n U. v.

(50)