各種平滑算法

1、算法：一、移動(dòng)窗口最小二乘多項(xiàng)式平滑(Savitzky-Golay Smoothing)假設(shè)數(shù)據(jù)(光譜或者是色譜等)為X,選定的平滑窗口大小為m(其必須為奇數(shù)，這里以7為例)，多項(xiàng)式次數(shù)為n,這里以3為例，當(dāng)前平滑的點(diǎn)為X0,前3個(gè)點(diǎn)分別記為：X-3,x_2,X-1 ,以及后三個(gè)點(diǎn)記為：Xi,X2,X30移動(dòng)窗口最小二乘多項(xiàng)式平滑就是利用中心點(diǎn)以及其前3個(gè)點(diǎn)和后3個(gè)點(diǎn)進(jìn)行最小二乘擬和。每一個(gè)點(diǎn)可以表示為不同的多項(xiàng)式的結(jié)果，從而7個(gè)點(diǎn)可以表示成為含有n+1(下面的例子是4個(gè))個(gè)未知數(shù)，m(例子中為7)個(gè)方程的方程組：X =bb(-3)b2*( -3)2b3*(-3)3= b-3bi9b2-27b

2、Xu =bb(-2)b2*( -2)2炳*(-2)3= b0-2tv4b2-8b3X- =bo+b!*( 1) + b2*( T)2+b3*(1)3=加bi+b2 b3X= b。bi*( -0) b2*(-0)2b3*(-0)3= bXi=bbi*(1)b2*(1)2b3*(1)3= bb2b3(1)X2 =b0bi*(2)b2*(2)2b32)3= b2tv 4b28b3X3= bbi*(3)b2*(3)2b3*(3)3= b3bi9b227b3對(duì)于上述方程的求解，采用最小二乘法。利用線性代數(shù)中的矩陣知識(shí), 為下面矩陣形式：即：A * b=x. . 一 .,_. 、 - . .一 、 -*因

3、而米用取小二乘法您算，礙到一個(gè)b b的解析解b :b*=( A危A )-1* At*x從而得到這個(gè)方程組的最小二乘解為:線性方程可以表示成1111111-3-2-1012941014x_iX。43 9 27X2X3Jlx3將求出來的b b*代入方程或者就可以求出平滑之后的數(shù)據(jù)點(diǎn)。實(shí)際上，如果將方程5求得的b b*代入方程1或者2之后得到如下7個(gè)方程：1 .Xw =一 *(39x-3+8x-2-4x-i-4x0+Xi+4x2-2x3)421 5x = *(8x-3+ 19X-2+ 16x-1+6x0-4xi-7x2+4x3)421.、Xj = 一*(-4x-3+ 16X-2+ 19X-1+12X

4、0+2X1-4X2+X3) 42X0 =*(-4x-3+6x-2+ 12x-1+ 14x0+ 12x1+6x2-4x3) 42x =*(x-3-4x-2+2x-1+12x0+19x1+ 16x2-4x3) 421 5X2 =*(4x-3-7x-2-4X-1+6X0+16X1+19x2+8x3) 421 X3 =一*(-2x-3+4x-2+X-1-4X0-4x1+8x2+39x3)42從這個(gè)里面我們可以發(fā)現(xiàn)，它們其實(shí)都是這個(gè)窗口內(nèi)部各個(gè)點(diǎn)的線性組合，即7個(gè)點(diǎn)由不同的權(quán)值進(jìn)行加權(quán)而得，對(duì)于我們需要的點(diǎn)X。也是由7個(gè)點(diǎn)加權(quán)而得。因此從本質(zhì)上說，移動(dòng)窗口多項(xiàng)式平滑其實(shí)就是利用窗口內(nèi)部各個(gè)點(diǎn)之間的加權(quán)來

5、計(jì)算平滑后的新值。計(jì)算過程中，中間局部我們只需要X0這個(gè)點(diǎn)的值即可，即從第四個(gè)點(diǎn)開始僅需要計(jì)算X0這個(gè)點(diǎn)的值。而對(duì)于開始的三個(gè)點(diǎn)和最后的三個(gè)點(diǎn)，沒有很好的處理方法，因此我們還是利用式子6來計(jì)算：開始的三個(gè)點(diǎn)用6式中的X-3 ,X-2,X-1計(jì)算式計(jì)算，最后的三個(gè)點(diǎn)用 式中的X1,X2,X3計(jì)算式計(jì)算。詳細(xì)解釋也可見分析化學(xué)手冊(cè)第十分冊(cè)。(-691821189-6 )15.5-16.75-14.5014.516.75-5.5一633.750-2.25-3-2.2503.751-1.751.751.750-1.75-1.751.75X0XiX2(5)(6)文二、粗糙懲罰(Roughness Pen

6、alty Smoothing)粗糙懲罰其實(shí)為了克服最小二乘法不穩(wěn)健而引入的一個(gè)方法。設(shè)平滑后的各個(gè)點(diǎn)為y*(i),最小二乘法的目標(biāo)函數(shù)是想讓最后的結(jié)果與原始數(shù)據(jù)之間的差異最小：2然而在實(shí)際情況中，如果有很多異常點(diǎn)的話， 這個(gè)標(biāo)準(zhǔn)并不能代表我們模型的準(zhǔn)確性，有時(shí)候反而會(huì)產(chǎn)生非常大的誤差，比方說色譜中如果噪聲水平很高的話，平滑效果并不好。因此,Silverman在1994出版的一本書中提出了粗糙懲罰算法，其就是在最小二乘目標(biāo)函數(shù)后面加上一個(gè)懲罰項(xiàng)：n2min (y(i) - y*(i)(；：2f (x)2d(x)i日式中，7、是懲罰系數(shù)，其越大，那么說明對(duì)這個(gè)數(shù)據(jù)點(diǎn)的懲罰越嚴(yán)重。后面的積分項(xiàng)是對(duì)函

7、數(shù)在x處的求二次導(dǎo)(這里的x并不是我們的數(shù)據(jù)點(diǎn)x(i),這個(gè)也就是高等數(shù)學(xué)里面的曲線的 曲率?，F(xiàn)在的問題是如何優(yōu)化這個(gè)目標(biāo)函數(shù)？目標(biāo)函數(shù)中前一個(gè)式子就是最小二乘擬和，可以通過回歸得到(同SG平滑)，而后面的積分式，由于f(x)很難得到。實(shí)際上，這個(gè)目標(biāo)函數(shù)是一個(gè)優(yōu)化問題，可將其轉(zhuǎn)化為線 性代數(shù)進(jìn)行求解。已經(jīng)證明了，如果函數(shù)f(x)可以通過立方樣條表示，那么可以通過一系列的變換得到如下的算式：(2f (x)2d(x) = yKy*(9)其中K通過下面的表達(dá)式求得:K = QRH(10)對(duì)于色譜或者光譜來講，由于是等間距采樣的，故可以得到Q和R的表達(dá)式如下:一2102-2(11)-2n*min (

8、y(i) - y(i)i =1(8)12/31/60000、1/62/31/600001/62/31/600R =00a1/6a2/39+ .0a09(12)00002/31/600001/62/3 j其中Q是一個(gè)n*(n-2)的一個(gè)矩陣，R是一個(gè)(n-2)*(n-2)的一個(gè)方形矩陣。利用上面兩個(gè)式子(11)和(12)代入方程(10)可以求出K,再代入方程(8)經(jīng)過變換之后，目 標(biāo)函數(shù)變?yōu)椋簄2 (y(i)-y*(i)(2f(x)2d(x)i V+ *+ *+ *=yt*y-2yt*y y*( I K)* y(13)求S的最小值。經(jīng)過變換可以發(fā)現(xiàn)，當(dāng):*_1y = (I K)1* y(14)的

9、時(shí)候，S可以取最小，這樣就求得了平滑函數(shù)的表達(dá)式。但是其中丸應(yīng)該如何判斷呢？在分析化學(xué)手冊(cè)第十分冊(cè)中提到了采用去一法交互檢驗(yàn)來選擇參數(shù)人，即：nCV = n、i=1*2(y(i) - y (%1-A()(15)其中A(舄)是矩陣A=(I+舄K)的第i個(gè)對(duì)角元素。通過代入不同的 人值可以得到不同 的CV值，在兀變化范圍之內(nèi)選擇CV值最小時(shí)對(duì)應(yīng)的值作為參數(shù)代入(14)式，就得到了平 滑后的函數(shù)。上面的兩個(gè)式子中，kern(x)就是kernel函數(shù)，是kernel平滑的核心，其有三種表達(dá)式, 分別代表不同的加權(quán)函數(shù)：10.5 x1均勻函數(shù)：kern(x)= 0 otherwise10. 7 5(1x

10、2) x苴1三、kernel平滑方法Kernel平滑方法在各種數(shù)據(jù)方法處理書中介紹得非常多，其本質(zhì)上和SG平滑一樣，采用加權(quán)函數(shù)，利用窗口內(nèi)部各個(gè)點(diǎn)的信息，來擬和當(dāng)前點(diǎn)，從而取得平滑效果。因此，kernel平滑方法的變化主要是變化在加權(quán)函數(shù)上面，而加權(quán)函數(shù)在kernel平滑以及書上介紹都是采用術(shù)語估計(jì)器(estimator)來表示。同樣假設(shè)平滑后的數(shù)據(jù)表示為y*(i)，帶寬為h,加權(quán)函數(shù)用S(t)表示，那么對(duì)于在某一時(shí)刻t (比方某個(gè)色譜的保存時(shí)間)的擬和值其計(jì)算式為：n* _ _y(t Sj(t)y(j)j =1(16)在1964年Nadaraya, Watson提出了加權(quán)函數(shù) 月(i)的計(jì)

11、算式:Sj(t)kern(L - t)/h)n、 kern(L - t)/h)j(17)后Gasser和Muller對(duì)其進(jìn)行了改良，下面提供的算法是改良后的算法，而Nadaraya-Watson加權(quán)函數(shù)可以直接代入式(16)進(jìn)行計(jì)算，因此不于詳細(xì)介紹。Gasser和Muller提供的加權(quán)函數(shù)表達(dá)式為：1Sj(t)=h,* t(j)t(j 1),U-tker n( ) h(18)二次函數(shù)：ken X =) S0otherwise高斯函數(shù)：ken x i(21/)2exup2(/ 2)利用不同的kernel函數(shù)會(huì)得到不同的平滑效果。在式子(17)中，t*=(tj+tj)/2 ,1j vn,t0=t

12、1,t；=tn。經(jīng)過變換后得到一個(gè)Gasser-Muller二次kernel加權(quán)函數(shù)為:1 33、1i3rj(t) - rj(t) -3rj(t) - rj(t)/ Sj(t) = 40otherwise式中，式(18)是通過式(17)和二次kernel用?？梢姡訖?quán)函數(shù)實(shí)質(zhì)上就是利用在帶寬加權(quán)，而對(duì)于在h范圍之外的點(diǎn)，不考慮其影響，即其權(quán)值為0,S是一個(gè)n*n，帶寬為h/2的帶狀矩陣，每一行代表一個(gè)點(diǎn)的加權(quán)信息。帶寬h的選擇一直以來都是數(shù)據(jù)處理學(xué)家們一個(gè)研究的課題，已經(jīng)提出非常多的方法來對(duì)其進(jìn)行自動(dòng)選擇，但是到現(xiàn)在還沒有一個(gè)可以解決所有問題的解決方案，因此對(duì)其算法也不予介紹。而且本人比擬信奉

13、的一句話是J.O. Ramsay和B.W. Silverman在他們的書functional data analysis中提到的h的選擇，任何一種算法其實(shí)還比不上人眼的判斷，如果我 們?cè)嚵藥讉€(gè)h值，覺得其中一個(gè)已經(jīng)到達(dá)了我們的要求了，那我們就無需花費(fèi)大量的心思在算法上面，把其變成完美(Our own view is that trying out a variety of values of h and inspectingthe consequences graphically remains a suitable means of resolving the bandwidth select