算法-Read.docVIP

算法： 一、移动窗口最小二乘多项式平滑(Savitzky-Golay Smoothing) 假设数据(光谱或者是色谱等)为x，选定的平滑窗口大小为m (其必须为奇数，这里以7为例)，多项式次数为n，这里以3为例，当前平滑的点为x0，前3个点分别记为：x-3，x-2，x-1，以及后三个点记为：x1，x2，x3。 移动窗口最小二乘多项式平滑就是利用中心点以及其前3个点和后3个点进行最小二乘拟和。每一个点可以表示为不同的多项式的结果，从而7个点可以表示成为含有n+1(下面的例子是4个)个未知数，m(例子中为7)个方程的方程组： (1) 对于上述方程的求解，采用最小二乘法。利用线性代数中的矩阵知识，线性方程可以表示成为下面矩阵形式： (2) 即： A*b=x (3) 因而采用最小二乘法运算，得到一个b的解析解 b*： b*=(At*A)-1*At*x (4) 从而得到这个方程组的最小二乘解为： (5) 将求出来的b*代入方程(1)或者(2)就可以求出平滑之后的数据点。实际上，如果将方程(5)求得的b*代入方程(1)或者(2)之后得到如下7个方程： (6) 从这个里面我们可以发现，它们其实都是这个窗口内部各个点的线性组合，即7个点由不同的权值进行加权而得，对于我们需要的点x0也是由7个点加权而得。因此从本质上说，移动窗口多项式平滑其实就是利用窗口内部各个点之间的加权来计算平滑后的新值。 计算过程中，中间部分我们只需要x0这个点的值即可，即从第四个点开始仅需要计算x0这个点的值。而对于开始的三个点和最后的三个点，没有很好的处理办法，因此我们还是利用式子(6)来计算：开始的三个点用(6)式中的x-3，x-2，x-1计算式计算，最后的三个点用(6)式中的x1，x2，x3计算式计算。 详细解释也可见分析化学手册第十分册。 二、粗糙惩罚 (Roughness Penalty Smoothing) 粗糙惩罚其实为了克服最小二乘法不稳健而引入的一个方法。设平滑后的各个点为y*(i)，最小二乘法的目标函数是想让最后的结果与原始数据之间的差别最小： (7) 然而在实际情况中，如果有很多异常点的话，这个标准并不能代表我们模型的准确性，有时候反而会产生非常大的误差，比如说色谱中如果噪声水平很高的话，平滑效果并不好。因此，Silverman在1994出版的一本书中提出了粗糙惩罚算法，其就是在最小二乘目标函数后面加上一个惩罚项： (8) 式中，是惩罚系数，其越大，则说明对这个数据点的惩罚越严重。后面的积分项是对函数在x处的求二次导(这里的x并不是我们的数据点x(i))，这个也就是高等数学里面的曲线的曲率。现在的问题是如何优化这个目标函数？ 目标函数中前一个式子就是最小二乘拟和，可以通过回归得到(同SG平滑)，而后面的积分式，由于很难得到。实际上，这个目标函数是一个优化问题，可将其转化为线性代数进行求解。已经证明了，如果函数f(x)可以通过立方样条表示，则可以通过一系列的变换得到如下的算式： (9) 其中K通过下面的表达式求得： (10) 对于色谱或者光谱来讲，由于是等间距采样的，故可以得到Q和R的表达式如下： (11) (12) 其中Q是一个n*(n-2)的一个矩阵，R是一个(n-2)*(n-2)的一个方形矩阵。 利用上面两个式子(11)和(12)代入方程(10)可以求出K，再代入方程(8)经过变换之后，目标函数变为： (13) 求S的最小值。经过变换可以发现，当： (14) 的时候，S可以取最小，这样就求得了平滑函数的表达式。 但是其中应该如何判断呢？在分析化学手册第十分册中提到了采用去一法交互检验来选择参数，即： (15) 其中是矩阵A=(I+K)的第i个对角元素。通过代入不同的值可以得到不同的CV值，在变化范围之内选择CV值最小时对应的值作为参数代入(14)式，就得到了平滑后的函数。 三、kernel平滑方法 Kernel平滑方法在各种数据方法处理书中介绍得非常多，其本质上和SG平滑一样，采用加权函