CN106096646A - 一种支持向量回归机模型选择方法 - Google Patents

Abstract

本发明公开了一种支持向量回归机模型选择方法，针对支持向量回归机的模型选择问题提出了一种新的基于混合核函数和容积卡尔曼滤波的支持向量回归机模型选择方法。选择混合核函数作为支持向量机的核函数，把混合核函数的组合系数嵌入到核函数参数与回归参数组成的超参数状态向量当中，从而将模型选择问题转换成一个非线性系统的状态估计问题，然后基于高性能的容积卡尔曼滤波进行超参数估计。仿真实验表明，本文所提方法与单个核函数的容积卡尔曼滤波支持向量回归机模型选择方法、遗传算法相比，该方法得到的决策回归函数具有更大的泛化能力，预测精度更高。

Description

一种支持向量回归机模型选择方法

技术领域

本发明涉及一种支持向量机相关技术，尤其涉及一种支持向量回归机模型选择方法。

背景技术

1995年Corinna Cortes和Vapnik等首先提出了支持向量机(Support VectorMachine，SVM)，它是建立在统计学习理论基础之上的通用学习方法。SVM由于其强大的非线性处理能力和泛化能力等特点，它在解决小样本、非线性及高维模式识别中表现出许多特有的优势，已经成功的解决了各种非线性和非可分机器学习问题，因此被广泛的应用在分类和回归问题中。但是，支持向量机方法也存在一些亟待解决和完善的问题，主要有：支持向量机分类理论是针对两类分类问题提出的，然而，现实世界的分类问题，如舰船识别、字体识别、人脸识别等都属于多分类的范畴，因此，如何用支持向量机更有效的解决多分类问题；支持向量机在二次寻优过程中要进行大量的矩阵运算，多数情况下，寻优算法占用了大部分的算法时间，这就使得存储空间和和计算时间成了求解二次规划问题的瓶颈，如何解决这一瓶颈问题；核函数以及核参数对支持向量机的分类拟合能力起着关键性作用，如何确定核函数以及最优的核参数以保证算法的有效性。

发明内容

本发明的目的就在于为了解决上述问题而提供一种支持向量回归机模型选择方法。

本发明通过以下技术方案来实现上述目的：

本发明包括参数滤波估计模型的建立、系统设计、预测输出函数、容积卡尔曼滤波、MKF-CKF-SVR模型参数选择；

参数滤波估计模型的建立：建立如下超参数非线性系统

γ(k)＝γ(k-1)+w(k) (1)

y(k)＝h(γ(k))+v(k) (2)

其中，γ(k)是超参数状态向量，y(k)是观测输出，过程噪声w(k)和观测噪声v(k)均是均值为零的高斯白噪声，且方差分别为Q和R；

由于待求的最优超参数可以看做是固定不变的，所以可以建立式(1)所示的关于超参数的线性状态方程，其次对于任何一个状态向量γ(k)，经过LIBSVM训练预测之后每一个原始数据都有一个预测输出，故可建立式(2)所示的非线性观测方程；为了CKF算法的运行，需要对系统模型加入人工过程白噪声和观测白噪声；

系统设计：支持向量回归机：

支持向量回归机的最终目的是找到一个回归函数f：R^D→R，使得

y＝f(x)＝w^Tφ(x)+b (3)

其中，φ(x)是一个将数据x从低维映射到高维特征空间的函数；w是一个权重向量，并且b是一个上下平移的数值；标准支持向量回归机采用ε-不灵敏函数，假设所有训练数据在精度ε下用线性函数拟合；这时，把问题转化为求优化目标函数最小化问题：

式中，ξ_i,是松弛因子，当拟合有误差时，ξ_i,都大于0，误差不存在时为0，优化函数第一项使拟合函数更为平坦，从而提高泛化能力；第二项为减小误差；常数C＞0表示对超出误差ε的样本的惩罚程度；

将核函数组合系数p1，p2、局部核函数的参数、全局核函数的参数以及惩罚参数C作为支持向量回归机的超参数γ，令k₁为局部核函数的核参数向量，k₂为全局核函数的核参数向量，则模型(1)中的超参数状态向量γ＝[p1,p2,k1,k2,C]^T；

求解(4)式的凸二次优化问题，引入Lagrange乘子α_i,支持向量回归机的原始问题(4)转换成如下对偶形式：

通过求解该对偶问题得到原始问题的解从而构造决策函数；将核函数K(x,x′)代替目标函数(1)中的內积(x_i·x_j)，则得到决策函数为：

其中，按如下方式计算：选择开区间中的或若选到的是则

若选到的是则

预测输出函数：假设支持向量回归机的原始样本数据集合为D＝{(x_i,y_i)|i∈I}，

其中指标集合I＝{1,2,...,N}，y_i为数据的目标向量，运用k-折交叉验证方法将样本数据分成k组，即

D_j＝{(x_i,y_i)|i∈I_j} (9)

其中j∈{1,2,...,k}，并且所有组的指标集I_j满足I₁∪I_j∪…∪I_k＝I，所有组的数据集D_j满足D₁∪D₂∪…∪D_k＝D；在每一次支持向量回归的迭代运算中，使用其中任意一组数据D_p用作预测，剩下的k-1组数据作为训练数据库，给定初始的超参数γ₀利用LIBSVM^[25]训练支持向量回归机；设此时训练结果为和则此时的决策函数为：

其中，

将数据组D_p代入式(11)，即可得到D_p的预测输出值

分别将数据组D_i，i∈{1,2,...,k}作为预测数据组，其余的数据组D₁,...,D_i-1,D_i+1,...,D_k作为支持向量回归机训练数据组，经过k-折交叉验证回归预测之后，样本数据集D中的每一个数据有且仅有一个预测输出值；故对于超参数向量γ，可定义如下预测输出函数：

y＝h(γ) (12)

其中，y＝(y(1),y(2),...,y(N))^T；

容积卡尔曼滤波：基于公式(1)、(2)所建立的超参数模型以及预测输出函数(12)，下面给出容积卡尔曼滤波算法部分的主要步骤，包括2个过程，即时间更新过程和测量更新过程：

时间更新：

1)假设在k时刻状态误差协方差矩阵已知，按下式分解

P(k-1|k-1)＝S(k-1|k-1)S^T(k-1|k-1) (13)

2)按下式计算容积点(i＝1,2,…,m)

$X_i(k-1|k-1)=S(k-1|k-1)\mathrm{\ξ}\left(i\right)+\widehat{\mathrm{\γ}}(k-1|k-1)---\left(14\right)$

其中，m＝2n_x，如若状态维数为2，则[1]∈R²，表示集合点且[1]_i表示该集合中的第i个向量点；

3)计算传播容积点(i＝1,2,…,m)

$X_i^{*}\left(k\right|k-1)=f\left(X_i\right(k|k-1\left)\right)---\left(15\right)$

4)计算一步状态预测

$\widehat{\mathrm{\γ}}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i^{*}\left(k\right|k-1)---\left(16\right)$

5)一步预测误差协方差阵为

$P\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i^{*}\left(k\right|k-1)X_i^{*T}\left(k\right|k-1)-\widehat{\mathrm{\γ}}\left(k\right|k-1){\widehat{\mathrm{\γ}}}^T\left(k\right|k-1)+Q\left(k\right)---\left(17\right)$

测量更新

1)按下式分解一步预测误差协方差阵

P(k|k-1)＝S(k|k-1)S^T(k|k-1) (18)

2)计算容积点(i＝1,2,…,m)

$X_i\left(k\right|k-1)=S\left(k\right|k-1)\mathrm{\ξ}\left(i\right)+\widehat{\mathrm{\γ}}\left(k\right|k-1)---\left(19\right)$

3)基于预测输出函数(14)-(17)计算传播容积点

Y_i(k|k-1)＝h(X_i(k|k-1)) (20)

4)一步测量预测为

$\hat{y}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{Y}_i\left(k\right|k-1)---\left(21\right)$

5)新息协方差矩阵为

$P_{yy}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{Y}_i\left(k\right|k-1)Y_i^T\left(k\right|k-1)-\hat{y}\left(k\right|k-1){\hat{y}}^T\left(k\right|k-1)+R\left(k\right)---\left(22\right)$

6)计算互协方差矩阵

$P_{xy}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i\left(k\right|k-1)Y_i^T\left(k\right|k-1)-\widehat{\mathrm{\γ}}\left(k\right|k-1){\hat{y}}^T\left(k\right|k-1)---\left(23\right)$

7)计算增益矩阵

$K\left(k\right)=P_{xy}\left(k\right|k-1)P_{yy}^{-1}\left(k\right|k-1)---\left(24\right)$

8)更新估计状态

$\widehat{\mathrm{\γ}}\left(k\right|k)=\widehat{\mathrm{\γ}}\left(k\right|k-1)+K\left(k\right)\left(y\right(k)-\hat{y}(k|k-1\left)\right)---\left(25\right)$

9)状态误差协方差矩阵为

P(k|k)＝P(k|k-1)-K(k)P_yy(k|k-1)K^T(k) (31)

Remark 2：由于CKF算法运用了径向积分和球形积分方法，该算法本身比UKF算法具有更高的估计精度；从容积卡尔曼滤波算法对超参数的估计过程来看，主要是在测量更新步骤(20)式中，需要将LIBSVM训练得到的预测输出函数嵌入到传播容积点的计算当中；

MKF-CKF-SVR模型参数选择

在超参数系统(1)-(2)中，观测向量y(k)的真实值在每一次迭代中都是不变的，且就是原始样本数据的目标值向量y(k)＝(y₁,y₂,...,y_N)^T，故可根据观测向量的真实值y(k)及预测输出值对超参数状态向量γ做最优状态估计，使得真实值与预测输出值的方差达到最小；MKF-CKF-SVR算法同样是包括两个过程，即时间更新过程和测量更新过程：

时间更新：

由于该更新过程是对状态的预测更新，且状态方程是线性已知的，故可根据容积卡尔曼滤波算法的时间跟新过程进行MKF-CKF-SVR算法的时间更新；

测量更新：

在测量更新的过程中需要用到预测输出函数，故不能直接利用CKF算法的公式(17)进行传播容积点的计算；需要由超参数状态向量γ(k)，利用LIBSVM训练数据集，再预测输出具体的MKF-CKF-SVR算法步骤：

初始化：

1.1获得原始数据集D，选择混合核函数模型，将组合系数嵌入到超参数状态向量γ当中，设置初始超参数状态值γ₀；

1.2运用k折交叉验证法首先将原始数据集分成k组D₁,D₂,…,D_k

While(超参数状态值不满足设定条件)do；

时间更新：

1.3利用公式(12)-(17)对超参数状态进行时间更新

测量更新：

1.4根据公式(17)分解一步预测误差协方差阵

1.5由用公式(18)计算容积点

1.6轮流将k组数据中的一组作为测试集，其它k-1组作为训练集，基于LIBSVM算法对数据集进行训练，最后预测输出

1.7基于预测利用公式(21)计算一步测量预测

1.8运用公式(21)-(25)继续进行测量更新

End while

End

Remark 3：基于混合核函数和容积卡尔曼滤波算法的SVR算法将混合核函数的组合系数与核参数及惩罚参数C作为超参数状态向量，然后利用k-折交叉验证法基于LIBSVM对数据集进行预测输出，最终用CKF算法迭代计算最优的的超参数状态向量；实际上，整个MKF-CKF-SVR算法的过程就是在迭代寻找最优的状态向量γ，使得样本真实目标值y(k)与支持向量回归机的预测输出之间的误差方差最小。

本发明的有益效果在于：

本发明是一种支持向量回归机模型选择方法，与现有技术相比，本发明将混合核函数的组合系数嵌入到核函数参数与回归参数组成的超参数状态向量当中，基于LIBSVM对原始数据集做预测输出，然后用容积卡尔曼滤波对超参数进行自动调整估计。最后以预测太阳黑子数平滑月均值为实验证明了基于本文所提方法得到的超参数能够使得支持向量回归机的泛化能力更强，预测精度更高。

附图说明

图1是本发明的MKF-CKF-SVR超参数调整结构图；

图2是本发明的MKF-CKF-SVR算法的回归预测图；

图3是本发明的RBF-GA-SVR算法的回归预测图。

具体实施方式

下面结合附图对本发明作进一步说明：

如图1所示：MKF-CKF-SVR模型参数选择

下面推导支持向量回归机模型超参数的选择方法，并给出所提算法的具体步骤。设计如图1所示的超参数调整系统，首先利用k-折交叉验证法将原始数据集分成k组，选择局部核函数和全局核函数确定混合核函数，基于该混合核函数用k个子LIBSVM训练此数据集，并将它的预测输出嵌入到容积卡尔曼滤波器当中，将模型的超参数作为系统的状态向量，则整个超参数的调整问题就可作为一个非线性动态系统的滤波估计问题。

本发明包括参数滤波估计模型的建立、系统设计、预测输出函数、容积卡尔曼滤波、MKF-CKF-SVR模型参数选择；

参数滤波估计模型的建立：建立如下超参数非线性系统

γ(k)＝γ(k-1)+w(k) (1)

y(k)＝h(γ(k))+v(k) (2)

其中，γ(k)是超参数状态向量，y(k)是观测输出，过程噪声w(k)和观测噪声v(k)均是均值为零的高斯白噪声，且方差分别为Q和R；

由于待求的最优超参数可以看做是固定不变的，所以可以建立式(1)所示的关于超参数的线性状态方程，其次对于任何一个状态向量γ(k)，经过LIBSVM训练预测之后每一个原始数据都有一个预测输出，故可建立式(2)所示的非线性观测方程；为了CKF算法的运行，需要对系统模型加入人工过程白噪声和观测白噪声；

系统设计：支持向量回归机：

支持向量回归机的最终目的是找到一个回归函数f：R^D→R，使得

y＝f(x)＝w^Tφ(x)+b (3)

其中，φ(x)是一个将数据x从低维映射到高维特征空间的函数；w是一个权重向量，并且b是一个上下平移的数值；标准支持向量回归机采用ε-不灵敏函数，假设所有训练数据在精度ε下用线性函数拟合；这时，把问题转化为求优化目标函数最小化问题：

式中，ξ_i,是松弛因子，当拟合有误差时，ξ_i,都大于0，误差不存在时为0，优化函数第一项使拟合函数更为平坦，从而提高泛化能力；第二项为减小误差；常数C＞0表示对超出误差ε的样本的惩罚程度；

将核函数组合系数p1，p2、局部核函数的参数、全局核函数的参数以及惩罚参数C作为支持向量回归机的超参数γ，令k₁为局部核函数的核参数向量，k₂为全局核函数的核参数向量，则模型(1)中的超参数状态向量γ＝[p1,p2,k1,k2,C]^T；

求解(4)式的凸二次优化问题，引入Lagrange乘子α_i,支持向量回归机的原始问题(4)转换成如下对偶形式：

通过求解该对偶问题得到原始问题的解从而构造决策函数；将核函数K(x,x′)代替目标函数(1)中的內积(x_i·x_j)，则得到决策函数为：

其中，按如下方式计算：选择开区间中的或若选到的是则

若选到的是则

预测输出函数：假设支持向量回归机的原始样本数据集合为D＝{(x_i,y_i)|i∈I}，

其中指标集合I＝{1,2,...,N}，y_i为数据的目标向量，运用k-折交叉验证方法将样本数据分成k组，即

D_j＝{(x_i,y_i)|i∈I_j} (9)

其中j∈{1,2,...,k}，并且所有组的指标集I_j满足I₁∪I_j∪…∪I_k＝I，所有组的数据集D_j满足D₁∪D₂∪…∪D_k＝D；在每一次支持向量回归的迭代运算中，使用其中任意一组数据D_p用作预测，剩下的k-1组数据作为训练数据库，给定初始的超参数γ₀利用LIBSVM^[25]训练支持向量回归机；设此时训练结果为和则此时的决策函数为：

其中，

将数据组D_p代入式(11)，即可得到D_p的预测输出值

分别将数据组D_i，i∈{1,2,...,k}作为预测数据组，其余的数据组D₁,...,D_i-1,D_i+1,...,D_k作为支持向量回归机训练数据组，经过k-折交叉验证回归预测之后，样本数据集D中的每一个数据有且仅有一个预测输出值；故对于超参数向量γ，可定义如下预测输出函数：

y＝h(γ) (12)

其中，y＝(y(1),y(2),...,y(N))^T；

容积卡尔曼滤波：基于公式(1)、(2)所建立的超参数模型以及预测输出函数(12)，下面给出容积卡尔曼滤波算法部分的主要步骤，包括2个过程，即时间更新过程和测量更新过程：

时间更新：

6)假设在k时刻状态误差协方差矩阵已知，按下式分解

P(k-1|k-1)＝S(k-1|k-1)S^T(k-1|k-1) (13)

7)按下式计算容积点(i＝1,2,…,m)

$X_i(k-1|k-1)=S(k-1|k-1)\mathrm{\ξ}\left(i\right)+\widehat{\mathrm{\γ}}(k-1|k-1)---\left(14\right)$

其中，m＝2n_x，如若状态维数为2，则[1]∈R²，表示集合点且[1]_i表示该集合中的第i个向量点；

8)计算传播容积点(i＝1,2,…,m)

$X_i^{*}\left(k\right|k-1)=f\left(X_i\right(k|k-1\left)\right)---\left(15\right)$

9)计算一步状态预测

$\widehat{\mathrm{\γ}}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i^{*}\left(k\right|k-1)---\left(16\right)$

10)一步预测误差协方差阵为

$P\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i^{*}\left(k\right|k-1)X_i^{*T}\left(k\right|k-1)-\widehat{\mathrm{\γ}}\left(k\right|k-1){\widehat{\mathrm{\γ}}}^T\left(k\right|k-1)+Q\left(k\right)---\left(17\right)$

测量更新

10)按下式分解一步预测误差协方差阵

P(k|k-1)＝S(k|k-1)S^T(k|k-1) (18)

11)计算容积点(i＝1,2,…,m)

$X_i\left(k\right|k-1)=S\left(k\right|k-1)\mathrm{\ξ}\left(i\right)+\widehat{\mathrm{\γ}}\left(k\right|k-1)---\left(19\right)$

12)基于预测输出函数(14)-(17)计算传播容积点

Y_i(k|k-1)＝h(X_i(k|k-1)) (20)

13)一步测量预测为

$\hat{y}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{Y}_i\left(k\right|k-1)---\left(21\right)$

14)新息协方差矩阵为

$P_{yy}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{Y}_i\left(k\right|k-1)Y_i^T\left(k\right|k-1)-\hat{y}\left(k\right|k-1){\hat{y}}^T\left(k\right|k-1)+R\left(k\right)---\left(22\right)$

15)计算互协方差矩阵

$P_{xy}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i\left(k\right|k-1)Y_i^T\left(k\right|k-1)-\widehat{\mathrm{\γ}}\left(k\right|k-1){\hat{y}}^T\left(k\right|k-1)---\left(23\right)$

16)计算增益矩阵

$K\left(k\right)=P_{xy}\left(k\right|k-1)P_{yy}^{-1}\left(k\right|k-1)---\left(24\right)$

17)更新估计状态

$\widehat{\mathrm{\γ}}\left(k\right|k)=\widehat{\mathrm{\γ}}\left(k\right|k-1)+K\left(k\right)\left(y\right(k)-\hat{y}(k|k-1\left)\right)---\left(25\right)$

18)状态误差协方差矩阵为

P(k|k)＝P(k|k-1)-K(k)P_yy(k|k-1)K^T(k) (31)

Remark 2：由于CKF算法运用了径向积分和球形积分方法，该算法本身比UKF算法具有更高的估计精度；从容积卡尔曼滤波算法对超参数的估计过程来看，主要是在测量更新步骤(20)式中，需要将LIBSVM训练得到的预测输出函数嵌入到传播容积点的计算当中；

MKF-CKF-SVR模型参数选择

在超参数系统(1)-(2)中，观测向量y(k)的真实值在每一次迭代中都是不变的，且就是原始样本数据的目标值向量y(k)＝(y₁,y₂,...,y_N)^T，故可根据观测向量的真实值y(k)及预测输出值对超参数状态向量γ做最优状态估计，使得真实值与预测输出值的方差达到最小；MKF-CKF-SVR算法同样是包括两个过程，即时间更新过程和测量更新过程：

时间更新：

由于该更新过程是对状态的预测更新，且状态方程是线性已知的，故可根据容积卡尔曼滤波算法的时间跟新过程进行MKF-CKF-SVR算法的时间更新；

测量更新：

在测量更新的过程中需要用到预测输出函数，故不能直接利用CKF算法的公式(17)进行传播容积点的计算；需要由超参数状态向量γ(k)，利用LIBSVM训练数据集，再预测输出具体的MKF-CKF-SVR算法步骤：

初始化：

1.1获得原始数据集D，选择混合核函数模型，将组合系数嵌入到超参数状态向量γ当中，设置初始超参数状态值γ₀；

1.2运用k折交叉验证法首先将原始数据集分成k组D₁,D₂,…,D_k

While(超参数状态值不满足设定条件)do；

时间更新：

1.3利用公式(12)-(17)对超参数状态进行时间更新

测量更新：

1.4根据公式(17)分解一步预测误差协方差阵

1.5由用公式(18)计算容积点

1.6轮流将k组数据中的一组作为测试集，其它k-1组作为训练集，基于LIBSVM算法对数据集进行训练，最后预测输出

1.7基于预测利用公式(21)计算一步测量预测

1.8运用公式(21)-(25)继续进行测量更新

End while

End

Remark 3：基于混合核函数和容积卡尔曼滤波算法的SVR算法将混合核函数的组合系数与核参数及惩罚参数C作为超参数状态向量，然后利用k-折交叉验证法基于LIBSVM对数据集进行预测输出，最终用CKF算法迭代计算最优的的超参数状态向量；实际上，整个MKF-CKF-SVR算法的过程就是在迭代寻找最优的状态向量γ，使得样本真实目标值y(k)与支持向量回归机的预测输出之间的误差方差最小。

算法分析

在基于局部核函数的的支持向量机的算法中，主要是相互之间距离较近的数据点对核函数的值产生影响，而在基于全局核函数的支持向量机算法中，主要是数据集中相互之间的距离较远的点对核函数的值产生影响。在求解实际问题时仅仅采用单一的核函数具有一定的局限性，往往不能使得得到的支持向量回归机性能最优。而混合核函数将两种不同类型的核函数组合之后使得新核函数既具有良好的学习能力又具有较好的推广能力。这就恰好体现了支持向量机在学习能力和推广能力之间两者兼顾的思想。但对于如何选择混合核函数的组合系数仍是个难题，通常采用的方法是专家的先验知识核基于交叉验证的方法，但这些方法都不一定能给出最优的选择参数。本文将构造核函数的组合参数p1，p2与局部核函数的参数、全局核函数的参数以及惩罚参数C放在一起作为支持向量回归机的超参数，事实上，我们可以完全将此时的超参数向量看作是混合核函数的核参数，然后基于高精度的容积卡尔曼滤波算法对超参数向量整体进行状态估计得到支持向量回归机最优的超参数值。

仿真例子

预测太阳黑子数平滑月均值的实验：将利用太阳活动1974年1月至2004年12月的太阳黑子数平滑月均值来预测太阳活动2005年1月至2014年11月的太阳黑子数平滑月均值。试验中的数据来自太阳影响数据分析中心。所有仿真均采用5-折交叉验证。本仿真混合核函数的局部核函数选择RBF核函数，全局核函数选择Sigmoid核函数。则MKF-CKF-SVR算法的超参数向量基于RBF核函数的遗传支持向量回归机算法(RBF-GA-SVR)的超参数向量γ＝[σ,C]^T。回归预测结果对比如图2和表2所示。

表2超参数估计结果表

从图2、图3可以看出，MKF-CKF-SVR算法与RBF-GA-SVR算法相比，能更好的对原始数据集进行拟合，预测精度更高。从表2的预测样本误差均方误差来看，MKF-CKF-SVR算法的预测样本误差均方误差大大小于RBF-GA-SVR算法，且基于MKF-CKF-SVR算法给出的局部核参数值σ远小于基于RBF-GA-SVR算法，因此，MKF-CKF-SVR算法泛化能力更强。

结论

本发明将混合核函数的组合系数嵌入到核函数参数与回归参数组成的超参数状态向量当中，基于LIBSVM对原始数据集做预测输出，然后用容积卡尔曼滤波对超参数进行自动调整估计。最后以预测太阳黑子数平滑月均值为实验证明了基于本文所提方法得到的超参数能够使得支持向量回归机的泛化能力更强，预测精度更高。

以上显示和描述了本发明的基本原理和主要特征及本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (1)

1.一种支持向量回归机模型选择方法，其特征在于：包括参数滤波估计模型的建立、系统设计、预测输出函数、容积卡尔曼滤波、MKF-CKF-SVR模型参数选择；

参数滤波估计模型的建立：建立如下超参数非线性系统

γ(k)＝γ(k-1)+w(k) (1)

y(k)＝h(γ(k))+v(k) (2)

其中，γ(k)是超参数状态向量，y(k)是观测输出，过程噪声w(k)和观测噪声v(k)均是均值为零的高斯白噪声，且方差分别为Q和R；

由于待求的最优超参数可以看做是固定不变的，所以可以建立式(1)所示的关于超参数的线性状态方程，其次对于任何一个状态向量γ(k)，经过LIBSVM训练预测之后每一个原始数据都有一个预测输出，故可建立式(2)所示的非线性观测方程；为了CKF算法的运行，需要对系统模型加入人工过程白噪声和观测白噪声；

系统设计：支持向量回归机：

支持向量回归机的最终目的是找到一个回归函数f：R^D→R，使得

y＝f(x)＝w^Tφ(x)+b (3)

其中，φ(x)是一个将数据x从低维映射到高维特征空间的函数；w是一个权重向量，并且b是一个上下平移的数值；标准支持向量回归机采用ε-不灵敏函数，假设所有训练数据在精度ε下用线性函数拟合；这时，把问题转化为求优化目标函数最小化问题：

式中，ξ_i,是松弛因子，当拟合有误差时，ξ_i,都大于0，误差不存在时为0，优化函数第一项使拟合函数更为平坦，从而提高泛化能力；第二项为减小误差；常数C＞0表示对超出误差ε的样本的惩罚程度；

将核函数组合系数p1，p2、局部核函数的参数、全局核函数的参数以及惩罚参数C作为支持向量回归机的超参数γ，令k₁为局部核函数的核参数向量，k₂为全局核函数的核参数向量，则模型(1)中的超参数状态向量γ＝[p1,p2,k1,k2,C]^T；

求解(4)式的凸二次优化问题，引入Lagrange乘子α_i,支持向量回归机的原始问题(4)转换成如下对偶形式：

通过求解该对偶问题得到原始问题的解从而构造决策函数；将核函数K(x,x')代替目标函数(1)中的內积(x_i·x_j)，则得到决策函数为：

其中，按如下方式计算：选择开区间中的或若选到的是则

若选到的是则

预测输出函数：假设支持向量回归机的原始样本数据集合为D＝{(x_i,y_i)|i∈I}，其中指标集合I＝{1,2,...,N}，y_i为数据的目标向量，运用k-折交叉验证方法将样本数据分成k组，即

D_j＝{(x_i,y_i)|i∈I_j} (9)

其中j∈{1,2,...,k}，并且所有组的指标集I_j满足I₁∪I_j∪…∪I_k＝I，所有组的数据集D_j满足D₁∪D₂∪…∪D_k＝D；在每一次支持向量回归的迭代运算中，使用其中任意一组数据D_p用作预测，剩下的k-1组数据作为训练数据库，给定初始的超参数γ₀利用LIBSVM^[25]训练支持向量回归机；设此时训练结果为和则此时的决策函数为：

其中，

将数据组D_p代入式(11)，即可得到D_p的预测输出值

分别将数据组D_i，i∈{1,2,...,k}作为预测数据组，其余的数据组D₁,...,D_i-1,D_i+1,...,D_k作为支持向量回归机训练数据组，经过k-折交叉验证回归预测之后，样本数据集D中的每一个数据有且仅有一个预测输出值；故对于超参数向量γ，可定义如下预测输出函数：

y＝h(γ) (12)

其中，y＝(y(1),y(2),...,y(N))^T；

容积卡尔曼滤波：基于公式(1)、(2)所建立的超参数模型以及预测输出函数(12)，下面给出容积卡尔曼滤波算法部分的主要步骤，包括2个过程，即时间更新过程和测量更新过程：

时间更新：

1)假设在k时刻状态误差协方差矩阵已知，按下式分解

P(k-1|k-1)＝S(k-1|k-1)S^T(k-1|k-1) (13)

2)按下式计算容积点(i＝1,2,…,m)

$X_i(k-1|k-1)=S(k-1|k-1)\mathrm{\ξ}\left(i\right)+\widehat{\mathrm{\γ}}(k-1|k-1)---\left(14\right)$

其中，m＝2n_x，如若状态维数为2，则[1]∈R²，表示集合点且[1]_i表示该集合中的第i个向量点；

3)计算传播容积点(i＝1,2,…,m)

$X_i^{*}\left(k\right|k-1)=f\left(X_i\right(k|k-1\left)\right)---\left(15\right)$

4)计算一步状态预测

$\widehat{\mathrm{\γ}}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i^{*}\left(k\right|k-1)---\left(16\right)$

5)一步预测误差协方差阵为

$P\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i^{*}\left(k\right|k-1)X_i^{*T}\left(k\right|k-1)-\widehat{\mathrm{\γ}}\left(k\right|k-1){\widehat{\mathrm{\γ}}}^T\left(k\right|k-1)+Q\left(k\right)---\left(17\right)$

测量更新

1)按下式分解一步预测误差协方差阵

P(k|k-1)＝S(k|k-1)S^T(k|k-1) (18)

2)计算容积点(i＝1,2,…,m)

$X_i\left(k\right|k-1)=S\left(k\right|k-1)\mathrm{\ξ}\left(i\right)+\widehat{\mathrm{\γ}}\left(k\right|k-1)---\left(19\right)$

3)基于预测输出函数(14)-(17)计算传播容积点

Y_i(k|k-1)＝h(X_i(k|k-1)) (20)

4)一步测量预测为

$\hat{y}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{Y}_i\left(k\right|k-1)---\left(21\right)$

5)新息协方差矩阵为

$P_{yy}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{Y}_i\left(k\right|k-1)Y_i^T\left(k\right|k-1)-\hat{y}\left(k\right|k-1){\hat{y}}^T\left(k\right|k-1)+R\left(k\right)---\left(22\right)$

6)计算互协方差矩阵

$P_{xy}\left(k\right|k-1)=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{X}_i\left(k\right|k-1)Y_i^T\left(k\right|k-1)-\widehat{\mathrm{\γ}}\left(k\right|k-1){\hat{y}}^T\left(k\right|k-1)---\left(23\right)$

7)计算增益矩阵

$K\left(k\right)=P_{xy}\left(k\right|k-1)P_{yy}^{-1}\left(k\right|k-1)---\left(24\right)$

8)更新估计状态

$\widehat{\mathrm{\γ}}\left(k\right|k)=\widehat{\mathrm{\γ}}\left(k\right|k-1)+K\left(k\right)\left(y\right(k)-\hat{y}(k|k-1\left)\right)---\left(25\right)$

9)状态误差协方差矩阵为

P(k|k)＝P(k|k-1)-K(k)P_yy(k|k-1)K^T(k) (31)

Remark 2：由于CKF算法运用了径向积分和球形积分方法，该算法本身比UKF算法具有更高的估计精度；从容积卡尔曼滤波算法对超参数的估计过程来看，主要是在测量更新步骤(20)式中，需要将LIBSVM训练得到的预测输出函数嵌入到传播容积点的计算当中；

MKF-CKF-SVR模型参数选择

在超参数系统(1)-(2)中，观测向量y(k)的真实值在每一次迭代中都是不变的，且就是原始样本数据的目标值向量y(k)＝(y₁,y₂,...,y_N)^T，故可根据观测向量的真实值y(k)及预测输出值对超参数状态向量γ做最优状态估计，使得真实值与预测输出值的方差达到最小；MKF-CKF-SVR算法同样是包括两个过程，即时间更新过程和测量更新过程：

时间更新：

由于该更新过程是对状态的预测更新，且状态方程是线性已知的，故可根据容积卡尔曼滤波算法的时间跟新过程进行MKF-CKF-SVR算法的时间更新；

测量更新：

在测量更新的过程中需要用到预测输出函数，故不能直接利用CKF算法的公式(17)进行传播容积点的计算；需要由超参数状态向量γ(k)，利用LIBSVM训练数据集，再预测输出具体的MKF-CKF-SVR算法步骤：

初始化：

1.1获得原始数据集D，选择混合核函数模型，将组合系数嵌入到超参数状态向量γ当中，设置初始超参数状态值γ₀；

1.2运用k折交叉验证法首先将原始数据集分成k组D₁,D₂,…,D_k

While(超参数状态值不满足设定条件)do；

时间更新：

1.3利用公式(12)-(17)对超参数状态进行时间更新

测量更新：

1.4根据公式(17)分解一步预测误差协方差阵

1.5由用公式(18)计算容积点

1.6轮流将k组数据中的一组作为测试集，其它k-1组作为训练集，基于LIBSVM算法对数据集进行训练，最后预测输出

1.7基于预测利用公式(21)计算一步测量预测

1.8运用公式(21)-(25)继续进行测量更新

End while

End

Remark3：基于混合核函数和容积卡尔曼滤波算法的SVR算法将混合核函数的组合系数与核参数及惩罚参数C作为超参数状态向量，然后利用k-折交叉验证法基于LIBSVM对数据集进行预测输出，最终用CKF算法迭代计算最优的的超参数状态向量；实际上，整个MKF-CKF-SVR算法的过程就是在迭代寻找最优的状态向量γ，使得样本真实目标值y(k)与支持向量回归机的预测输出之间的误差方差最小。