各种平滑算法

算法：一、移动窗口最小二乘多项式平滑(Savitzky-GolaySmoothing)假设数据(光谱或者是色谱等)为x，选定的平滑窗口大小为m(其必须为奇数，这里以7为例)，多项式次数为n，这里以3为例，当前平滑的点为x0，前3个点分别记为：x-3，x-2，x-1，以及后三个点记为：x1，x2，x3。移动窗口最小二乘多项式平滑就是利用中心点以及其前3个点和后3个点进行最小二乘拟和。每一个点可以表示为不同的多项式的结果，从而7个点可以表示成为含有n+1(下面的例子是4个)个未知数，m(例子中为7)个方程的方程组：23301230123232012301232310123012323001230210123*(3)*(3)*(3)3927*(2)*(2)*(2)248*(1)*(1)*(1)*(0)*(0)*(0)*(1)*(1)*(1)xbbbbbbbbxbbbbbbbbxbbbbbbbbxbbbbbxbbbb301232320123012323301230123*(2)*(2)*(2)248*(3)*(3)*(3)3927bbbbxbbbbbbbbxbbbbbbbb(1)对于上述方程的求解，采用最小二乘法。利用线性代数中的矩阵知识，线性方程可以表示成为下面矩阵形式：320110213231392712481111*10001111124813927xxbxbxbxbxx(2)即：A*b=x(3)因而采用最小二乘法运算，得到一个b的解析解b*：b*=(At*A)-1*At*x(4)从而得到这个方程组的最小二乘解为：32011021323-691821189-65.5-16.75-14.5014.516.75-5.51*3.750-2.25-3-2.2503.7563-1.751.751.750-1.75-1.751.75xxbxbxbxbxx(5)将求出来的b*代入方程(1)或者(2)就可以求出平滑之后的数据点。实际上，如果将方程(5)求得的b*代入方程(1)或者(2)之后得到如下7个方程：3-3-2-101232-3-2-101231-3-2-101230-3-2-101231-3-2-1011*(39x+8x-4x-4x+x+4x-2x)421*(8x+19x+16x+6x-4x-7x+4x)421*(-4x+16x+19x+12x+2x-4x+x)421*(-4x+6x+12x+14x+12x+6x-4x)421*(x-4x+2x+12x+19x+42xxxxx232-3-2-101233-3-2-1012316x-4x)1*(4x-7x-4x+6x+16x+19x+8x)421*(-2x+4x+x-4x-4x+8x+39x)42xx(6)从这个里面我们可以发现，它们其实都是这个窗口内部各个点的线性组合，即7个点由不同的权值进行加权而得，对于我们需要的点x0也是由7个点加权而得。因此从本质上说，移动窗口多项式平滑其实就是利用窗口内部各个点之间的加权来计算平滑后的新值。计算过程中，中间部分我们只需要x0这个点的值即可，即从第四个点开始仅需要计算x0这个点的值。而对于开始的三个点和最后的三个点，没有很好的处理办法，因此我们还是利用式子(6)来计算：开始的三个点用(6)式中的x-3，x-2，x-1计算式计算，最后的三个点用(6)式中的x1，x2，x3计算式计算。详细解释也可见分析化学手册第十分册。二、粗糙惩罚(RoughnessPenaltySmoothing)粗糙惩罚其实为了克服最小二乘法不稳健而引入的一个方法。设平滑后的各个点为y*(i)，最小二乘法的目标函数是想让最后的结果与原始数据之间的差别最小：2*1min(()())niyiyi(7)然而在实际情况中，如果有很多异常点的话，这个标准并不能代表我们模型的准确性，有时候反而会产生非常大的误差，比如说色谱中如果噪声水平很高的话，平滑效果并不好。因此，Silverman在1994出版的一本书中提出了粗糙惩罚算法，其就是在最小二乘目标函数后面加上一个惩罚项：2*221min(()())(())()niyiyifxdx(8)式中，是惩罚系数，其越大，则说明对这个数据点的惩罚越严重。后面的积分项是对函数在x处的求二次导(这里的x并不是我们的数据点x(i))，这个也就是高等数学里面的曲线的曲率。现在的问题是如何优化这个目标函数？目标函数中前一个式子就是最小二乘拟和，可以通过回归得到(同SG平滑)，而后面的积分式，由于2()fx很难得到。实际上，这个目标函数是一个优化问题，可将其转化为线性代数进行求解。已经证明了，如果函数f(x)可以通过立方样条表示，则可以通过一系列的变换得到如下的算式：22**(())()tfxdxyKy(9)其中K通过下面的表达式求得：1tKQRQ(10)对于色谱或者光谱来讲，由于是等间距采样的，故可以得到Q和R的表达式如下：10000002100000121000001210000012000000001210000001200000001Q(11)2/31/600001/62/31/600001/62/31/600001/62/30000002/31/600001/62/3R(12)其中Q是一个n*(n-2)的一个矩阵，R是一个(n-2)*(n-2)的一个方形矩阵。利用上面两个式子(11)和(12)代入方程(10)可以求出K，再代入方程(8)经过变换之后，目标函数变为：2*221***(()())(())()*2**()*nitttSyiyifxdxyyyyyIKy(13)求S的最小值。经过变换可以发现，当：*1()*yIKy(14)的时候，S可以取最小，这样就求得了平滑函数的表达式。但是其中应该如何判断呢？在分析化学手册第十分册中提到了采用去一法交互检验来选择参数，即：2*11()()()1()niiyiyiCVnA(15)其中()iA是矩阵A=(I+K)的第i个对角元素。通过代入不同的值可以得到不同的CV值，在变化范围之内选择CV值最小时对应的值作为参数代入(14)式，就得到了平滑后的函数。三、kernel平滑方法Kernel平滑方法在各种数据方法处理书中介绍得非常多，其本质上和SG平滑一样，采用加权函数，利用窗口内部各个点的信息，来拟和当前点，从而取得平滑效果。因此，kernel平滑方法的变化主要是变化在加权函数上面，而加权函数在kernel平滑以及书上介绍都是采用术语估计器(estimator)来表示。同样假设平滑后的数据表示为y*(i)，带宽为h，加权函数用S(t)表示，则对于在某一时刻t(比如某个色谱的保留时间)的拟和值其计算式为：*1()()()njjytStyj(16)在1964年Nadaraya，Watson提出了加权函数()jSi的计算式：1ker(()/)()ker(()/)jjnjjntthStntth(17)后Gasser和Muller对其进行了改进，下面提供的算法是改进后的算法，而Nadaraya-Watson加权函数可以直接代入式(16)进行计算，因此不于详细介绍。Gasser和Muller提供的加权函数表达式为：*()*(1)1()ker()jjtjtutStnhh(18)上面的两个式子中，kern(x)就是kernel函数，是kernel平滑的核心，其有三种表达式，分别代表不同的加权函数：均匀函数：0.51ker()0xnxotherwise二次函数：20.75(1)1ker()0xxnxotherwise高斯函数：1/22ker()(2)exp(/2)nxu利用不同的kernel函数会得到不同的平滑效果。在式子(17)中，*1()/2jjjttt，1j＜n，*01tt，*nntt。经过变换后得到一个Gasser-Muller二次kernel加权函数为：33111[3()()][3()()]()40jjjjjjrtrtrtrtrthStotherwise(19)式中，*()jjttrth式(18)是通过式(17)和二次kernel函数推出来的一个表达式，对上式我们可以直接利用。可见，加权函数实质上就是利用在带宽h(相似于移动窗口大小)内部的各个点的信息来加权，而对于在h范围之外的点，不考虑其影响，即其权值为0，从而得到的一个加权矩阵S是一个n*n，带宽为h/2的带状矩阵，每一行代表一个点的加权信息。带宽h的选择一直以来都是数据处理学家们一个研究的课题，已经提出非常多的方法来对其进行自动选择，但是到现在还没有一个可以解决所有问题的解决方案，因此对其算法也不予介绍。而且本人比较信奉的一句话是J.O.Ramsay和B.W.Silverman在他们的书functionaldataanalysis中提到的h的选择，任何一种算法其实还比不上人眼的判断，如果我们试了几个h值，觉得其中一个已经达到了我们的要求了，那我们就无需花费大量的心思在算法上面，把其变成完美(Ourownviewisthattryingoutavarietyofvaluesofhandinspectingtheconsequencesgraphicallyremainsasuitablemeansofresolvingthebandwidthselectionproblemformostpracticalproblems.)。对于kernel函数的详细介绍可见J.O.Ramsay和B.W.Silverman著的functionaldataanalysis，springer2005年出版。