分段多項式及樣條


本文永久鏈接:分段多項式及樣條 | ESL CN

寫在前面

最早接觸樣條是在科學計算的課程上,當時主要講了三次樣條及由此推出的三斜率方程組。

后來便在閱讀ESL(The Elements of Statistical Learning)這本書的時候,接觸了樣條,ESL中的介紹就更加系統,以及更加偏統計了。

第三次正式接觸樣條是金融風險管理的課程中,課程最后,張老師簡單介紹了樣條估計。

本篇博客主要是ESL的第五章第2節的翻譯,順帶提一句,我現在一直在維護ESL-CN這個網站,這是個ESL的翻譯項目,目前已完成書中大部分章節的翻譯,以及一小部分的代碼實現,歡迎大家訪問並提出寶貴的建議😄

分段多項式

直到5.7節我們都假設 \(X\) 為一維向量。分段多項式函數\(f(X)\)通過將\(X\)的定義域分成連續的區間,然后在每個區間內用單獨的多項式來表示\(f\)。圖5.1顯示了兩個簡單的分段多項式。第一個是分段常值,含有三個基函數:

\[h_1(X)=I(X<\xi_1),\;h_2(X)=I(\xi_1\le X\le \xi_2),\;h_3(X)=I(\xi_2\le X) \]

因為這些在不連續區域為正值,模型\(f(X)=\sum_{m=1}^3\beta_mh_m(X)\)等價於\(\hat\beta_m=\overline{Y}_m\),即為\(Y\)在第\(m\)個區域的均值。

圖5.1

圖5.1. 上面兩張圖顯示了對一些擬合數據的分段常值函數擬合。垂直虛線表示兩個結點\(\xi_1\)\(\xi_2\)。藍色曲線表示真正的函數,數據是通過函數加上高斯噪聲產生的。下兩張圖顯示了對同樣數據的分段線性函數擬合——上面板的圖沒有限制,而下面板的圖限制為在結點處連續——右上圖沒有限制,而左下圖限制為在結點處連續。右下圖顯示了分段線性的基函數,\(h_3(X)=(X-\xi_1)_+\),它在\(\xi_1\)處連續。黑色點表示樣本取值\(h_3(x_i),i=1,2,\ldots,N\)

右上圖顯示了分段線性擬合。需要三個額外的基函數:\(h_{m+3}=h_m(X)X,m=1,2,3\).除了特殊的情形,我們一般更想要第三個面板的圖,也是分段線性,但在間隔點上連續。這些連續性的約束導致在參數上的線性約束;舉個例子,\(f(\xi_1^-)=f(\xi_1^+)\)意味着\(\beta_1+\xi_1\beta_4=\beta_2+\xi_1\beta_5\).在這種情形下,因為存在兩個約束,我們期減少了兩個參數,最終得到4個自由參數。

這種情形下更直接的方式是將約束結合起來的基函數:

\[h_1(X)=1,\;h_2(X)=X,\;h_3(X)=(X-\xi_1)_+,\;h_4(X)=(X-\xi_2)_+ \]

其中\(t_+\)記為正的部分。函數\(h_3\)的圖象如圖5.1顯示。我們經常偏好光滑函數,這個可以通過增加局部多項式的階數實現。圖5.2顯示了一系列對同樣數據進行的分段3次多項式擬合,增加了結點處的連續性的階數。右下圖的函數是連續的,且在結點處的一階微分和二階微分均連續。這稱為三次樣條。再加一階的連續性可以導致全局三次多項式。不難證明(練習5.1)下面的基函數表示結點為\(\xi_1\)\(\xi_2\)的三次樣條。

\[\begin{align} h_1(X)=1,\;h_3(X)=X^2,\;&h_5(X)=(X-\xi_1)_+^3\\ h_2(X)=X,\;h_4(X)=X^3,\;&h_6(X)=(X-\xi_2)_+^3 \end{align} \qquad (5.3) \]

這里6個基函數對應6維函數的線性空間。快速地確定參數個數:(3個區域)\(\times\) (每個區域4個參數)-(2個結點)\(\times\)(每個結點3個限制)=6.

weiya注
單獨考慮每個區域,則需要確定的參數有四個,\(\sum_{i=1}^4\beta_mh_m\),每個結點處,需要保持連續、一階微分連續、二階微分連續,所以每個結點3個限制。
從另一個角度看,因\(\sum_{i=1}^6\beta_mh_m\)需要確定的就是\(\beta_i\)這6個參數。

圖5.2

圖5.2. 一系列分段3次多項式擬合,增加了連續性的階數。

更一般地,結點為\(\xi_j,j=1,\ldots,K\)的order為\(M\)的樣條是order為\(M\)的分段多項式,而且有連續的\(M-2\)次微分。三次樣條的\(M=4\)。實際上圖5.1的分段常數函數是order為1的樣條,而連續的分段線性函數是order為2的樣條。同樣地,截斷(truncated-power)基的集合是

\[\begin{align} h_j(X)&=X^{j-1},j=1,\ldots,M\\ h_{M+\ell}(X)&=(X-\xi_\ell)_+^{M-1},\ell=1,\ldots,K \end{align} \]

據說三次樣條是人眼看不出結點不連續的最低階樣條。很少有更好的理由去選擇更高次的樣條,除非對光滑的微分感興趣。實際中,用得最多的order還是\(M=1,2,4\)

這些固定結點的樣條也稱作回歸樣條(regression splines)。我們需要選擇樣條的階數,結點的個數以及它們的位置。一種簡單方式是用基函數或自由度來參量化樣條族,並用觀測\(x_i\)來確定結點的位置。舉個例子,R語言中的命令bs(x,df=7)產生在x\(N\)個觀測點取值的三次樣條基函數,其中有\(7-3=4\)個內結點,內結點在x的(20,40,60和80)分位數處。(含四個結點的三次樣條是8個維度的。bs()函數默認忽略基函數里面的常數項,因為這樣的項一般包含在模型的其它項里面。)然而,也可以更明確地指出,`bs(x, degree=1, knots=c(0.2,0.4,0.6))``產生有三個內結點的線性樣條的基,並且返回一個\(N\times 4\)的矩陣。

因為特定階數以及結點序列的樣條函數的空間是向量空間,所以表示它們會有許多等價的基底(就像普通多項式一樣。)盡管truncated power基在概念上很簡單,但是數值計算時不是很吸引人:大的power會導致非常嚴重的舍入問題。在本章附錄中描述的B樣條的基即使在結點數\(K\)很大時也有很高的計算效率。

自然三次樣條

我們知道對數據的多項式擬合的行為在邊界處有不穩定的趨勢,而且外推法會很危險。樣條進一步惡化了這些問題。邊界結點之外的多項式擬合的表現比該區域對應的全局多項式擬合更野蠻。從最小二乘擬合的樣條函數的逐點方差(pairwise variance),可以很方便地看到這一點(更多細節見下一節的計算這些方差的例子)。圖5.3比較了不同模型的逐點方差。在邊界處的方差爆炸是顯而易見的,對於三次樣條這必要會更糟糕。

自然三次樣條(natural cubic spline)添加額外的限制,具體地,令邊界結點之外的函數是線性的。這樣減少了4個自由度(兩個邊界區域分別兩個限制條件),這四個自由度可以通過在內部區域取更多的結點花費掉。圖5.3用方差表示了這種權衡。在邊界附近需要在偏差上付出代價,但是假設邊界附近(不管怎樣,我們的信息很少)為線性函數通常是合理的考慮。

\(K\)個結點的自然三次樣條用\(K\)個基函數來表示。可以從三次樣條的出發,通過強加上邊界限制導出降維的基。舉個例子,從5.2節描述的truncated power序列基出發,我們得到(練習5.4):

\[N_1(X)=1,\;N_2(X)=X,\; N_{k+2}(X)=d_k(X)-d_{K-1}(X),\qquad (5.4) \]

其中,

\[d_k(X)=\frac{(X-\xi_k)_+^3-(X-\xi_K)_+^3}{\xi_K-\xi_k}\qquad (5.5) \]

可以看到當\(X\ge \xi_K\)時每個基函數的二階微分和三階微分均為0.

例子:南非心臟病

在4.4.2節我們對南非心臟病數據進行了線性邏輯斯蒂回歸擬合。這里我們以采用自然樣條的函數來探索非線性。模型的函數有如下形式:

\[\mathrm{logit}[Pr(chd\mid X)]=\theta_0+h_1(X_1)^T\theta_1+h_2(X_2)^T\theta_2+\cdots+h_p(X_p)^T\theta_p\qquad (5.6) \]

其中每個\(\theta_i\)都是乘以對應自然樣條基函數\(h_j\)的系數向量。

我們在模型中對每一項采用4個自然樣條基。舉個例子,\(X_1\)代表sbp,\(h_1(X_1)\)是包含四個基函數的基。這實際上表明有三個而非兩個內結點(在sbp的均勻分位數結點處取值),另外在數據端點有兩個邊界結點,因為我們把\(h_j\)的常數項提取出來了。

weiya注
下面說明對於自然三次樣條而言,基函數個數即為結點個數。
設有\(K\)個結點,則有\((K+1)\)個區域,每個區域參數為4個,每個內結點減掉3個參數,每個邊界點減掉兩個參數,則還剩下

$$
(K+1)\cdot 4-3K-2\times 2 = K
$$
也就是$K$個基函數。

因為 famhist是含兩個水平的因子,所以用一個二進制變量或者虛擬變量來編碼,而且它與擬合的模型中的單系數有關。

更簡潔地,我們將\(p\)維基函數(以及常數項)向量結合成一個向量\(h(X)\),則模型簡化為\(h(X)^T\theta\),總參數個數為\(df=1+\sum_{j=1}^pdf_j\),是每個組分中參數個數的和。每個基函數在\(N\)個樣本中分別取值,得到一個\(N\times df\)的基矩陣\(\mathbf H\).在這點上看,模型類似於其他的線性邏輯斯蒂回歸模型,應用的算法在4.4.1節描述。

我們采用向后逐步刪除過程,從模型中刪除項並且保持每個項的整體結構,而不是每次刪除一個系數。AIC統計量(7.5節)用來刪除項,並且在最后模型中保留下來的所有項如果被刪掉都會導致AIC增大(見表5.1).圖5.4顯示了通過逐步回歸選擇出的最終模型的圖象。對於每個變量\(X_j\),畫出的函數是\(\hat f_j(X_j)=h_j(X_j)^T\hat\theta_j\)。協方差矩陣\(Cov(\hat\theta)=\mathbf\Sigma\)通過\(\mathbf{\hat\Sigma=(H^TWH)^{-1}}\)來估計,其中\(\mathbf W\)為邏輯斯蒂回歸的對角元素構成的權重矩陣。因此\(v_j(X_j)=Var[\hat f_j(X_j)]=h_j(X_j)^T\mathbf{\hat \Sigma}\_{jj}h_j(X_j)\)\(\hat f_j\)的逐點方差函數,其中\(Cov(\hat\theta_j)=\hat{\mathbf\Sigma}\_{jj}\)\(\hat{\mathbf\Sigma}\)對應的子矩陣。圖5.4中每張圖的陰影區域由\(\hat f_j(X_j)\pm2\sqrt{v_j(X_j)}\)定義。

AIC統計量比似然比檢驗(偏差檢驗)更“寬容”(generous)。sbpobesity都被包含進模型中,而這兩個量都不在線性模型中。圖象解釋了為什么它們的貢獻本質上是非線性。這些影響乍看或許很奇怪,但這因為是回顧性數據的本質。這些指標有時是當病人患上心臟病后測出來的,而且在很多情形下他們已經受益於健康飲食和生活狀態,因此在obesitysbp值較低時會有明顯的增長。表5.1總結了部分模型的效果。

例子:音素識別

在這個例子中我們采用樣條來降低靈活性而非增大靈活性;這個應用是屬於一般的函數型(functional)建模。圖5.5的上圖顯示了在256個頻率下分別測量兩個音素“aa”和“ao”的15個對數周期圖。目標是應用這些數據對口語音素進行分類。選擇這兩個音素是因為它們很難被分開。

圖5.5 上圖顯示了對數周期圖作為15個例子頻率的函數,例子中每個音素“aa”和“ao”從總共695個“aa”和1022個“ao”中選取。每個對數周期圖在256個均勻的空間頻率處測量。下圖顯示了對數據進行極大似然擬合邏輯斯蒂回歸得到的系數(作為頻率的函數)的圖象,將256個對數周期圖值作為輸入。系數被限制為在紅色曲線中是光滑的,而在鋸齒狀灰色曲線中不受限制。

weiya注
周期圖(Periodogram):在信號處理中,周期圖是信號譜密度的估計。

輸入特征是長度為256的向量\(x\),我們可以看成是在頻率\(f\)的節點上取值的函數值\(X(f)\)向量。實際上,存在一個連續的類似信號,它是頻率的函數,而且我們有它的一個取樣版本。

圖5.5的下圖顯示了對從695個“aa”和1022個“ao”選出的1000個訓練樣本進行極大似然擬合得到的線性邏輯斯蒂回歸模型的系數。也作出了系數關於頻率的函數圖象,而且實際上我們可以根據下面的連續形式來思考模型

\[\mathrm{log}\frac{Pr(aa\mid X)}{Pr(ao\mid X)}=\int X(f)\beta(f)df\qquad (5.7) \]

可以用下式來近似

\[\sum\limits_{j=1}^{256}X(f_j)\beta(f_j)=\sum\limits_{j=1}^{256}x_j\beta_j\qquad (5.8) \]

系數計算出對比度函數(contrast functional),而且將會在頻域內有顯著的值,其中對數周期圖會區分這兩個類。

灰色曲線十分粗糙。因為輸入信號有相當強的正自相關性,這導致系數中的負自相關。另外,樣本大小僅僅為每個系數僅僅提供了4個有效的觀測。

類似這樣的應用允許自然的正則化。我們強制系數作為頻率的函數均勻變化。圖5.5中下圖的紅色曲線顯示了對這些數據應用這樣一個光滑參數曲線。我們看到低頻率的差異性很明顯。這個光滑不僅允許對相反的進行簡單的解讀,而且得到更加精確的分類器:

紅色光滑曲線可以應用非常簡單的自然三次樣條得到。我們可以將系數函數表達成樣條\(\beta(f)=\sum_{m=1}^Mh_m(f)\theta_m\)的展開。實際中這意味着\(\beta=\mathbf H\theta\),其中,\(\mathbf H\)\(p\times M\)三次樣條的基矩陣,定義在頻率集上。這里我們采用\(M=12\)個基函數,其中結點均勻分布在表示頻率的整數\(1,2,\ldots,256\)上。因為\(x^T\beta=x^T\mathbf H\theta\),我們可以簡單地將輸入特征\(x\)替換成濾波(filtered)形式\(x^*=\mathbf{H}^Tx\),並在\(x^*\)上通過線性邏輯斯蒂回歸擬合\(\theta\)。因此紅色曲線是\(\beta(f)=h(f)^T\hat\theta\)

友情鏈接:

  1. IR^p 中結構化局部回歸
  2. 徑向基函數和核
  3. 7.7 貝葉斯方法和BIC
  4. 12.3 支持向量機和核


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM