本发明涉及地理信息科学领域，尤其涉及一种基于AIC的逐步特征变量选择方法。
背景技术：
地理加权回归是近年来提出的一种新的空间分析方法，它是将空间结构嵌入线性回归模型中，能够有效处理回归分析中空间非平稳性现象。在地理空间分析当中，利用在n个不同的地理位置上获取的样本数据进行回归分析时，回归参数在很多情况下是随着地理位置变化的，为了反映回归参数的空间变化特征，因而需要采用地理加权回归的分析方法。
在地理加权回归(Geographically Weighted Regression，GWR)分析时，需要从多个特征变量中选择合适的变量，作为自变量与因变量一起建立模型。特征变量的选择对模型性能影响很大。在普通线性回归(Ordinarily Linearity Regression，OLR)中，常用的特征变量选择方法包括基于F检验的向前引入法、向后剔除法和逐步回归法。由于GWR在OLR的基础上考虑了各影响因素在空间非平稳性，上述方法无法直接应用于GWR。在GWR特征变量选择方法上，赵阳阳提出一种基于贪心算法的地理加权回归特征变量选择方法该方法基于Akaike信息法则(Akaike Information Criterion，AIC)，采用贪心算法筛选变量，提高了GWR中特征变量选择的精度。但是，由于贪心算法选择过程是单向的，一旦变量选入模型或删除后，无法再判断，得到的结果可能只是局部最优，而不是全局最优。
技术实现要素：
鉴于上述现有技术中存在的以上问题，本发明结合AIC准则和逐步筛选原理，提出一种基于AIC的逐步特征变量选择方法，解决GWR分析中无法建立全局最优模型的问题。
本发明提供的基于AIC的逐步特征变量选择方法包括以下步骤：
步骤1：设全部特征变量集为L，记作L＝{x1，x2，…，xd}，其中d大于等于2；随机选择L中的1个特征变量放入已选特征变量集S，令备选特征变量集U为子集S在L中的绝对补集，记作U＝{x1，x2，…，xk}，全局最小AIC值为AICw＝+∞，即初始赋值为无穷大；
步骤2：计算备选特征变量集U中所包含的备选特征变量的个数，记作k；判断k是否为0；若是，则执行步骤11；若否，则执行步骤3；
步骤3：对备选特征变量集U中的每个备选特征变量，作为新变量分别与已选特征变量集S、因变量y建立GWR模型，得到模型GWRu1，GWRu2，...，GWRuk；计算模型GWRu1，GWRu2，...，GWRuk当中每个模型的AIC值，获得{AICGWRu1，AICGWRu2，...，AICGWRuk}，其中AICGWRuk表示模型GWRuk的AIC值；
步骤4：求出{AICGWRu1，AICGWRu2，...，AICGWRuk}当中的最小值，记为当前最小AIC值AICc；
步骤5：判断步骤4中所述AICc是否比全局最小AIC值AICw小3；若是，执行步骤6；若否，随机从U中选择一个变量，删除该变量，执行步骤7；
步骤6：将计算出步骤4中所述AICc的GWR模型所对应的备选特征变量加入已选特征变量集S，将加入后的已选特征变量集S记作S＝{x1，x2，…，xm}；并从备选特征变量集U中删除该备选特征变量，将所述AICc的值赋值给AICw；然后执行步骤7；
步骤7：对已选特征变量集S中的m个变量，每次去除其中一个变量，利用S中的其它变量与因变量y组合建立GWR模型，得到模型GWRs1，GWRs2，...，GWRsm；其中GWRsm表示将S中的变量xm去除后，利用其它变量与因变量y组合而建立的GWR模型；并计算模型GWRs1，GWRs2，...，GWRsm当中每个模型的AIC值，获得{AICGWRs1，AICGWRs2，...，AICGWRsm}，其中AICGWRsm表示模型GWRsm的AIC值；
步骤8：求{AICGWRs1，AICGWRs2，...，AICGWRsm}当中的最小值，记为当前最小AIC值AICc；
步骤9：判断步骤8中所述AICc是否比全局最小AIC值AICw小3；若是，执行步骤10；若否，执行步骤2；
步骤10：将计算出步骤8中所述AICc的GWR模型对应的特征变量从已选特征变量集S中删除，将步骤8中所述AICc赋值给AICw，返回执行步骤2；
步骤11：循环结束，此时的已选特征变量集作为最优特征变量组合。
其中，优选的是，步骤3和步骤7中因变量y表示n个数据观测点的观测值，公式为：
其中，优选的是，步骤3中的GWRu1模型为：
GWRu2模型为：
GWRuk模型为：
其中xi1，xi2，……，xim表示数据点(ui，vi)处已选特征变量集S中的m个变量，xiu1，xiu2，……，xiuk表示数据点(ui，vi)处备选特征变量集U中的k个变量；βip(ui，vi)，(p＝0，1，…，m)是第i个数据点(ui，vi)处变量xip对应的未知参数，βiu1(ui，vi)、βiu2(ui，vi)……βiuk(ui，vi)是第i个数据点(ui，vi)处分别与变量xiu1，xiu2，……，xiuk对应的未知参数，ui表示第i个数据点的横坐标，vi表示第i个数据点的纵坐标；(ε1，ε2，…εn)为独立同分布的误差项，通常假定均值为零，方差为σ2；i的取值范围为1至n，表示所述n个数据观测点。
其中，优选的是，步骤7中的GWRs1模型为：
yi＝βi0(ui，vi)+βi2(ui，vi)xi2+βi3(ui，vi)xi3+…+βim(ui，vi)xim+εi (5)
GWRs2模型为：
yi＝βi0(ui，vi)+βi1(ui，vi)xi1+βi3(ui，vi)xi3+…+βim(ui，vi)xim+εi (6)
GWRsm模型为：
yi＝βi0(ui，vi)+βi1(ui，vi)xi1+βi2(ui，vi)xi2+…+βim-1(ui，vi)xi(m-1)+εi (7)
其中xi1，xi2，……，xi(m-1)，xim表示第i个数据点(ui，vi)处已选特征变量集S中的m个变量；
βi1(ui，vi)，βi2(ui，vi)，…，βim(ui，vi)是第i个数据点(ui，vi)处变量xi1，xi2，…，xim对应的未知参数，ui表示第i个数据点的横坐标，vi表示第i个数据点的纵坐标；(ε1，ε2，…εn)为独立同分布的误差项，通常假定均值为零，方差为σ2；i的取值范围为1至n，表示所述n个数据观测点。
其中，优选的是，步骤3和步骤7中每个模型的AIC值计算公式如下：
其中，n为数据观测点个数，tr(S)为地理加权回归帽子矩阵S的迹，是误差项估计的标准离差：
为n个数据观测点的估计值：
其中帽子矩阵S为：
X表示X的转置，X表示GWR模型中所有特征变量的n个数据点的观测值：公式(2)中的X为：
公式(3)中的X为：
公式(4)中的X为：
公式(5)中的X为：
公式(6)中的X为：
公式(7)中的X为：
Wi为权重矩阵：
其中wij为数据观测点i与数据观测点j的权重；采用高斯函数法计算权重：
wij＝exp(-(dij/b)2) (19)
dij表示点(ui，vi)与点(ui，vi)之间的距离；采用欧式距离计算：
b表示最优带宽，采用CV法求最优带宽：
表示在回归参数估计时不包含该点本身；设置b取值范围为(0，max(dij))，当CV最小时，b对应的值即为最优带宽。
步骤4中的AICc为：
AICc＝min(AICGWRu1，AICGWRu2，...，AICGWRuk) (22)
步骤8中的AICc为：
AICc＝min(AICGWRs1，AICGWRs2，...，AICGWRsm) (23)
与现有技术相比，本发明所述的方法定义了在GWR建模中特征变量的选择方法，采用逐步选择的方式，解决了单向方法易陷入局部最优的问题，能从全局角度选择最优组合，建立精确模型。
说明书附图
图1是本发明所述的一种基于AIC的逐步特征变量选择方法的流程图。
具体实施方式
下面通过实施例，对本发明的技术方案做进一步具体的说明。
本发明结合AIC准则和逐步筛选原理，提出一种基于AIC的逐步特征变量选择方法，解决GWR分析中无法建立全局最优模型的问题。
图1是本发明提供的基于AIC的逐步特征变量选择方法的流程图，可见，本发明具体包括以下步骤：
步骤1：设全部特征变量集为L，记作L＝{x1，x2，…，xd}其中d大于等于2；随机选择L中的1个特征变量放入已选特征变量集S，令备选特征变量集U为子集S在L中的绝对补集，记作U＝{x1，x2，…，xk}，全局最小AIC值为AICw＝+∞，即初始赋值为无穷大；
步骤2：计算备选特征变量集U中所包含的备选特征变量的个数，记作k；判断k是否为0；若是，则执行步骤11；若否，则执行步骤3；
步骤3：对备选特征变量集U中的每个备选特征变量，作为新变量分别与已选特征变量集S、因变量y建立GWR模型，得到模型GWRu1，GWRu2，...，GWRuk；计算模型GWRu1，GWRu2，...，GWRuk当中每个模型的AIC值，获得{AICGWRu1，AICGWRu2，...，AICGWRuk}，其中AICGWRuk表示模型GWRuk的AIC值；
步骤4：求出{AICGWRu1，AICGWRu2，...，AICGWRuk}当中的最小值，记为当前最小AIC值AICc；
步骤5：判断步骤4中所述AICc是否比全局最小AIC值AICw小3；若是，执行步骤6；若否，随机从U中选择一个变量，删除该变量，执行步骤7；
步骤6：将计算出步骤4中所述AICc的GWR模型所对应的备选特征变量加入已选特征变量集S，将加入后的已选特征变量集S记作S＝{x1，x2，…，xm}；并从备选特征变量集U中删除该备选特征变量，将所述AICc的值赋值给AICw；然后执行步骤7；
步骤7：对已选特征变量集S中的m个变量，每次去除其中一个变量，利用S中的其它变量与因变量y组合建立GWR模型，得到模型GWRs1，GWRs2，...，GWRsm，；其中GWRsm表示将S中的变量xm去除后，利用其它变量与因变量y组合而建立的GWR模型；并计算模型GWRs1，GWRs2，...，GWRsm当中每个模型的AIC值，获得{AICGWRs1，AICGWRs2，...，AICGWRsm}，其中AICGWRsm表示模型GWRsm的AIC值；
步骤8：求{AICGWRs1，AICGWRs2，...，AICGWRsm}当中的最小值，记为当前最小AIC值AICc；
步骤9：判断步骤8中所述AICc是否比全局最小AIC值AICw小3；若是，执行步骤10；若否，执行步骤2；
步骤10：将计算出步骤8中所述AICc的GWR模型对应的特征变量从已选特征变量集S中删除，将步骤8中所述AICc赋值给AICw，返回执行步骤2；
步骤11：循环结束，此时的已选特征变量集作为最优特征变量组合。
其中，步骤3和步骤7中因变量y表示n个数据观测点的观测值，公式为：
步骤3中的GWRu1模型为：
GWRu2模型为：
GWRuk模型为：
其中xi1，xi2，……，xim表示数据点(ui，vi)处已选特征变量集S中的m个变量，xiu1xiu2，……，xiuk表示数据点(ui，vi)处备选特征变量集U中的k个变量；βip(ui，vi)，(p＝0，1，…，m)是第i个数据点(ui，vi)处变量xip对应的未知参数，βiu1(ui，vi)、βiu2(ui，vi)……βiuk(ui，vi)是第i个数据点(ui，vi)处分别与变量xiu1，xiu2……，xiuk对应的未知参数，ui表示第i个数据点的横坐标，vi表示第i个数据点的纵坐标；(ε1，ε2，…εn)为独立同分布的误差项，通常假定均值为零，方差为σ2；i的取值范围为1至n，表示所述n个数据观测点。
其中，优选的是，步骤7中的GWRs1模型为：
yi＝βi0(ui，vi)+βi2(ui，vi)xi2+βi3(ui，vi)xi3+…+βim(ui，vi)xim+εi (5)
GWRs2模型为：
yi＝βi0(ui，vi)+βi1(ui，vi)xi1+βi3(ui，vi)xi3+…+βim(ui，vi)xim+εi (6)
GWRsm模型为：
yi＝βi0(ui，vi)+βi1(ui，vi)xi1+βi2(ui，vi)xi2+…+βim-1(ui，vi)xi(m-1)+εi (7)
其中xi1，xi2，……，xi(m-1)，xim表示第i个数据点(ui，vi)处已选特征变量集S中的m个变量；
βi1(ui，vi)，βi2(ui，vi)，…，βim(ui，vi)是第i个数据点(ui，vi)处变量xi1，xi2，…，xim对应的未知参数，ui表示第i个数据点的横坐标，vi表示第i个数据点的纵坐标；(ε1，ε2，…εn)为独立同分布的误差项，通常假定均值为零，方差为σ2；i的取值范围为1至n，表示所述n个数据观测点。
其中，优选的是，步骤3和步骤7中每个模型的AIC值计算公式如下：
其中，n为数据观测点个数，tr(S)为地理加权回归帽子矩阵S的迹，是误差项估计的标准离差：
为n个数据观测点的估计值：
其中帽子矩阵S为：
X表示X的转置，X表示GWR模型中所有特征变量的n个数据点的观测值：公式(2)中的X为：
公式(3)中的X为：
公式(4)中的X为：
公式(5)中的X为：
公式(6)中的X为：
公式(7)中的X为：
Wi为权重矩阵：
其中wij为数据观测点i与数据观测点j的权重；本专利采用高斯函数法计算权重：
wij＝exp(-(dij/b)2) (19)
dij表示点(ui，vi)与点(ui，vi)之间的距离；本专利采用欧式距离计算：
b表示最优带宽，本专利采用CV法求最优带宽，
表示在回归参数估计时不包含该点本身；设置b取值范围为(0，max(dij))，当CV最小时，b对应的值即为最优带宽。
步骤4中的AICc为：
AICc＝min(AICGWRu1，AICGWRu2，...，AICGWRuk) (22)
步骤8中的AICc为：
AICc＝min(AICGWRs1，AICGWRs2，...，AICGWRsm) (23)
与现有技术相比，本发明所述的方法定义了在GWR建模中特征变量的选择方法，采用逐步选择的方式，解决了单向方法易陷入局部最优的问题，能从全局角度选择最优组合，建立精确模型。
以上实施例仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。