部分线性模型的高斯径向基函数估计

Gaussian Radial Basis Function Estimation of Partially Linear Model

摘要: 部分线性模型是一种常用的现代统计模型，其同时具备参数与非参数回归的优点。我们基于高斯径向基函数估计部分线性模型的非线性部分，并给出估计过程中的超参数选择方法。在模拟仿真与实证分析中将高斯径向基函数与B样条进行对比，发现高斯径向基函数在部分线性模型中可以成为B样条的一种替代方法。

Abstract: The partially linear model is a commonly used modern statistical model with the advantages of both parametric and nonparametric regression. We estimate the nonlinear part of the partially linear model based on the Gaussian radial basis function, and give the hyperparameter selection method in the estimation process. The Gaussian radial basis function is compared with the B-spline in simulation and empirical analysis, and it is found that the Gaussian radial basis function can be an alternative method to the B-spline in the partially linear model.

1. 引言

部分线性模型是一种常用的统计模型，其将参数模型与非参数模型的优势相结合，用于解决具有复杂关系的数据分析问题。在部分线性模型中，解释变量可以包含线性部分和非线性部分，通过将线性和非线性部分相结合，部分线性模型能够更好地拟合数据并提高预测准确性。部分线性模型形式如下：

$Y_i={X^{\prime }}_i\beta +g\left(Z_i\right)+{ϵ}_i,i=1,2,\cdots ,n$

其中 $Y_i$ 是响应变量， $X_i$ 是线性部分解释变量， $\beta$ 是线性部分系数。 $Z_i$ 是非线性部分解释变量， $g$ 为未知函数。 ${ϵ}_i$ 为误差项， ${ϵ}_i$ 与 $X_i$ ， $Z_i$ 独立。令 $ϵ=\left({ϵ}_1,\cdots ,{ϵ}_n\right)\text{'}$ ， $ϵ$ 满足 $Var\left(ϵ\right)={\rho }^2{I}_n$ ， $I_n$ 为n阶单位阵。

由于部分线性模型性质良好，其相关应用与理论受到了众多关注。实际应用方面，1986年，Engle等 [1] 首次将部分线性模型应用于用电需求相关问题。1988年，Speckman [2] 运用该模型探究了漱口水的效果。1999年，Schmalensee和Stoker [3] 运用该模型研究了美国家庭汽油消费情况。理论性质方面，1986年，Heckman [4] 基于光滑样条估计部分线性模型非参数部分。1988年，Robinson [5] 研究了该模型非参数部分的核估计方法，并给出估计方法的大样本性质。2006年，Ma等 [6] 研究了异方差下该模型的加权估计方法，并证明了相合性。2021年，Liu和Yin [7] 研究了时间序列数据下部分线性模型的样条估计方法。Zhong等 [8] 基于深度神经网络估计部分线性Cox模型中的非线性部分，并给出了估计方法的大样本性质。Rodríguez等 [9] 研究了具有单调性约束部分线性模型的稳健估计方法，并通过蒙特卡洛随机模拟验证了方法的有效性。

经典的估计方法主要使用核与样条。我们考虑了一种新的基函数：高斯径向基函数。由于高斯径向基函数的良好性质，有许多学者对其进行了研究。2008年，Ando等 [10] 使用信息准则选取了正则化径向基函数神经网络的超参数。2014年，Lei等 [11] 利用局部Rademacher复杂度研究了径向基函数网络的泛化性能，获得了新的估计误差边界。2021年，Krzyżak和Niemann [12] 研究了归一化径向基函数网络的收敛性与收敛速率，讨论了径向基函数和网络参数的选择。2022年，Sosa与Buitrago [13] 运用径向基函数逼近时变系数模型的系数函数，并给出频率方法与贝叶斯估计方法。

本文的组织如下：第二章介绍基于高斯径向基函数的估计方法。第三章进行模拟仿真，并与B样条进行对比。第四章，将估计方法应用于实际数据集。第五章，对全文总结。

2. 估计方法

2.1. 模型估计方法

使用高斯径向基函数估计部分线性模型就是使用一组高斯径向基函数的线性组合去逼近未知函数g。一组高斯径向基函数形式如下：

$\exp \left(-\frac{{\left(z-c_1\right)}^2}{2{\sigma }^2}\right),\cdots ,\exp \left(-\frac{{\left(z-c_k\right)}^2}{2{\sigma }^2}\right),$

其中 ${\left\{c_j\right\}}_{j=1}^k$ 为高斯径向基函数的中心参数， $\sigma$ 为形状参数。由此有

$g(z)\approx {\displaystyle \underset{j=1}{\overset{k}{\sum }}{\gamma }_j}\exp \left(-\frac{{\left(z-c_j\right)}^2}{2{\sigma }^2}\right).$

则有

$Y_i\approx {X^{\prime }}_i\beta +{\displaystyle \underset{j=1}{\overset{k}{\sum }}{\gamma }_j}\exp \left(-\frac{{\left(Z_i-c_j\right)}^2}{2{\sigma }^2}\right)+{ϵ}_i,i=1,2,\cdots ,n,$

通过最小二乘求解参数，即

$\underset{\beta ,\gamma ,\sigma }{\min }{\left(Y_i-{X^{\prime }}_i\beta -{\displaystyle \underset{j=1}{\overset{k}{\sum }}{\gamma }_j}\exp \left(-\frac{{\left(Z_i-c_j\right)}^2}{2{\sigma }^2}\right)\right)}^2$ (1)

其中 $\gamma =\left({\gamma }_1,\cdots ,{\gamma }_k\right)$ 。以上的优化问题是一个非凸优化。记 $\mathbb{X}$ 为(1)的设计阵， $\mathbb{Y}={\left(Y_1,\cdots ,Y_t\right)}^{\prime }$ 。

Figure 1. Diagram of non-sparse point

图1. 数据点均匀分布图

2.2. 模型参数选择

目标函数(2)为一个非凸优化问题，无法直接应用优化算法求解。我们将该问题分三步求解，最终将转化为凸优化问题。

首先进行中心参数选择。中心参数选择分为两个部分：中心位置与中心个数。中心参数位置的选择方法主要有两类。第一类不依赖于 ${\left\{Z_i\right\}}_{i=1}^n$ 的分布，如等间距点与分位数点。这种方法是多项式样条的常用取点方式，高斯径向基函数在一维下也是适用的。这类取点方法有着很多优点，如计算成本低，理论性质好等。等间距点适用的 ${\left\{Z_i\right\}}_{i=1}^n$ 如图1，在z轴均匀分布。

第二类方法依赖于 ${\left\{Z_i\right\}}_{i=1}^n$ 的分布，如Kmeans聚类以及在 ${\left\{Z_i\right\}}_{i=1}^n$ 中随机抽取。第二类方法显然比第一类更加适合在样本点稀疏分布时取点。但其有着缺陷：随着中心参数的增加，中心参数之间的距离变化不易描述，甚至Kmeans聚类无法保证中心间的最小距离趋于0，且Kmeans聚类计算量过大。等间距点适用的 ${\left\{Z_i\right\}}_{i=1}^n$ 如图2，在z轴非均匀分布。本文主要考虑一维问题，所以采用等间距点。

Figure 2. Diagram of sparse point

图2. 数据点非均匀分布图

中心个数过多会产生过拟合，反之则欠拟合。样条类的方法根据大样本性质可以直接使用 $n^{1/5}$ 作为节点个数。高斯径向基函数无渐近理论参考，我们采用数据驱动的中心个数选择方法。中心个数选择通过如下信息准则

$\begin{array}{l}\text{AIC}=\log \left(\frac{\text{RSS}}{n}\right)+\frac{2K}{n};\\ \text{AICc}=\log \left(\frac{\text{RSS}}{n}\right)+\frac{2K}{n}+\frac{2\left(K+1\right)\left(K+2\right)}{n\left(n-K-2\right)};\\ \text{BIC}=\log \left(\frac{\text{RSS}}{n}\right)+\log \left(n\right)\frac{K}{n},\end{array}$

其中n为样本量，K为(1)待估参数个数，RSS为(1)的残差平方和最小值。使用了AIC，AICc，BIC等方法进行参数选择，其数值结果显示AIC是一种简单且效果良好的方法，因此采用AIC准则。信息准则是根据样本内的误差进行模型选择。若要使用样本外误差进行模型选择，需使用交叉验证以及广义交叉验证，可以参考 [14] 。

其次，位置参数选择。高斯径向基函数的逼近效果依赖于形状参数σ。选择方法可主要分为三类，第一类取决于节点的位置，可参考 [15] 。第二类由数据驱动选择，如AIC与交叉验证等。第三类同时依赖于位置与个数，可参考 [16] 。

Figure 3. Diagram of Gaussian radial basis functions under different σ

图3. 不同σ下的高斯径向基函数

以上三类方法有一种共同思想，当中心点稀疏时，σ应越大，从而使得高斯径向基函数趋于平稳，影响范围变大，基函数不同σ下的变化趋势如图3。结合以上思想，本文采取 $\sigma =M_0/\sqrt{k},M_0>0$ 。对于 $M_0$ 的选择可以考虑AIC，BIC，交叉验证等方法。本文为减少计算量直接采用 $M_0=1$ 。

最后求解问题(1)的系数。在给定中心参数与未知参数后，问题(1)为普通最小二乘问题，即有显示解。

3. 模拟仿真

本节中，我们对比基于高斯径向基函数与基于B样条的估计方法。考虑如下模型：

$Y_i=X_i\beta +g\left(Z_i\right)+{ϵ}_i,i=1,.2,\cdots ,n$ (2)

$X_i$ 与 $Z_i$ 服从 $U\left(0,1\right)$ ， ${ϵ}_i$ 服从 $N\left(0,0.25\right),\beta =2$ 。 $g\left(z\right)$ 考虑如下5种函数： $x,x^2,\sin \left(x\right),\exp \left(-x^2\right)$ 以及 $\log \left(1+x\right)$ ，将其分别简记为 $f_1,\cdots ,f_5$ 。

按照 $X_i,Z_i,{ϵ}_i$ 的分布进行独立抽样，样本量n为100，200，300，400，500。再利用模型(2)计算 $Y_i$ ，从而得到观测数据 $\left(X_i,Z_i,Y_i\right)$ 。基于该观测数据分别使用高斯径向函数与B样条估计未知函数 $g\left(z\right)$ 。取[0, 1]上等间距点 ${\left\{Xtes{t}_j\right\}}_{j=1}^{100},{\left\{Ztes{t}_j\right\}}_{j=1}^{100}$ ，基于 ${\left(Xtes{t}_j,Ztes{t}_j\right)}_{j=1}^{100}$ 与模型(2)计算 ${\left\{Ytes{t}_j\right\}}_{j=1}^{100}$ ，但计算过程中不加入扰动。使用 ${\left(Xtes{t}_j,Ztes{t}_j,Ytes{t}_j\right)}_{j=1}^{100}$ 与MAE比较两种方法，MAE公式如下：

$\text{MAE}=\frac{{\displaystyle {\sum }_{j=1}^{100}\left|\hat{Y}tes{t}_j-Ytes{t}_j\right|}}{100}$

其中 $\hat{Y}tes{t}_j$ 为估计方法的预测值。

Table 1. Mean of 500 MAEs of Gaussian radial basis functions

表1. 高斯径向基函数500次MAE均值

Table 2. Mean of 500 MAEs of Bspline

表2. B样条500次MAE均值

Figure 4. Box plot of 500 MAEs of $f_1$

图4. $f_1$ 500次MAE的箱线图

Figure 5. Box plot of 500 MAEs of $f_2,\cdots ,f_5$

图5. $f_2,\cdots ,f_5$ 500次MAE的箱线图

B样条节点与高斯径向基函数中心都采用等间距点，个数范围为2~10，使用AIC准则选取节点个数。为减少模拟随机性，重复模拟500次。使用500次试验MAE的均值衡量估计效果，试验结果如表1与表2。我们发现高斯径向基函数在5个函数上都优于B样条，且在 $f_4$ 上差距表现最为明显。同时基于高斯径向基函数的估计方法随着样本量的增大，MAE逐渐减小。这说明估计方法具有一致性。图4与图5为5个函数下500次MAE的箱线图，高斯径向基函数在5个函数上MAE的稳定性都高于B样条，这说明未知函数估计上高斯径向基函数更为稳定。高斯径向基函数在 $f_1,\cdots ,f_5$ 上表现也有所不同， $f_1,f_2$ 的异常值更多更大， $f_3,f_4,f_5$ 异常值相比则较少， $f_5$ 在样本量少时稳定性最高。

4. 实证研究

实证研究使用波士顿房价数据集，该数据集可以从R包MASS中获取。刘志伟和夏志明 [17] 也使用该数据集研究了半线性模型。我们基于数据集中medv，rm，lstat，ptratio，dis等变量进行数据分析，其中medv为被解释变量，其余为解释变量。medv为房价。

各变量间的散点图如图6。rm为每间住宅的平均房间数，容易发现rm与medv的散点图有线性的趋势，且呈现房间数上升房价上升的规律，这符合常理。lstat为低收入群比例，lsta与medv的散点图依然具有线性的趋势，且呈现低收入群比例上升房价降低的规律，数据是符合现实的。ptratio为城镇中的教师学生比例，ptratio与medv的散点图任然呈现线性关系，但随着ptratio变化的上升下降趋势并不显著。dis为距离5个波士顿的就业中心的加权距离，dis与medv的散点图既具有线性趋势，但又可能具有更复杂的函数关系。同时具有距离就业中心越远，房价越高的倾向。

由于dis与medv不一定呈现线性关系，计算相关系数用于判断变量间线性相关程度。相关系数表如表3。容易发现medv与rm，lstat，ptratio，相关系数较大，结合散点图判断medv与这些变量之间呈线性关系。medv与dis之间相关系数为0.250，数值较小，我们结合散点图推断medv与dis之间为非线性关系。

Figure 6. Scatterplot between variables

图6. 变量间散点图

Table 3. Table of correlation coefficients

表3. 相关系数表

根据变量间的线性与非线性关系，我们考虑建立如下部分线性模型：

$med{v}_i={\beta }_{0}r{m}_i+{\beta }_{1}lsta{t}_i+{\beta }_{2}ptrati{o}_i+g\left(di{s}_i\right)+{ϵ}_i,i=1,2,\cdots ,506.$

使用该模型对比高斯径向基函数与B样条，将数据10折，计算MAE用于对比两种方法的数值效果。高斯径向基函数与B样条的10折MAE如表4，均值分别为0.086与0.087，高斯径向基函数的表现优于B样条。这反应了高斯径向基函数在实证有着不错的效果，可以在实际数据中成为B样条的替代方法。

Table 4. 10 fold MAE

表4. 10折MAE

5. 结论

我们使用高斯径向基函数估计了部分线性模型中的非参数部分，分析了高斯径向基函数的超参数选择方法。同时，在5种不同的函数下对比了高斯径向基函数与B样条的数值效果，在波士顿房价数据集下对比了两种方法的实际表现。我们发现高斯径向基函数都优于B样条，由此认为高斯径向基函数可以成为B样条在部分线性模型中的一种替代方法。

参考文献