非参数回归模型样条估计量的分布

Distribution of Spline Estimators for Nonparametric Regression Models

摘要: 为探究非参数回归模型中非参数函数估计量的分布，本文在标准正态误差情形下首先得到了均值函数样条估计量的正态分布，然后得到了方差函数基于残差的样条估计量的渐近分布，并采用单个卡方变量线性函数来近似方差函数估计量的渐近分布。通过数值模拟验证了均值函数估计量的分布和方差函数估计量的渐近分布。

Abstract: To explore the distribution of nonparametric function estimator in the nonparametric regression model, we first obtain the normal distribution of the spline estimator of the mean function in the standard normal error casein this paper. Then the asymptotic distribution of the spline estimator of the variance function based on the residuals is obtained. And the linear function of individual chi-square variable is used to approximate the asymptotic distribution of the variance function es-timator. The distribution of the mean function estimator and the asymptotic distribution of the variance function estimator are illustrated by numerical simulations.

1. 引言

考虑如下非参数回归模型

$y_t=m\left(x_t\right)+\sigma \left(x_t\right){\epsilon }_t,t=1,2,\cdots ,n$， (1)

其中 $x_t=t/n$ 为固定设计点， $y_t$ 为被解释变量或响应变量， $m(\cdot )$ 与 $\sigma (\cdot )$ 分别为未知的均值函数和标准差函数， $\left\{{\epsilon }_t\right\}$ 为独立同分布的随机误差序列，且 $E{\epsilon }_t=0$， $E{\epsilon }_t^2=1$， $t=1,2,\cdots ,n$。

近些年来，国内外学者对非参数回归模型做了大量的研究 [1] - [7]。渐近分布作为估计量大样本性质中的一个重要方面，也引起了一些学者的兴趣。秦永松(1991) [8] 基于加权核估计得到了均值函数的估计量，并研究了其导数的渐近分布；Liang与Jing (2004) [9] 采用加权核估计的方法，基于负相关序列研究了非参数回归模型中未知函数估计量的逐点一致收敛性与渐近正态性；Jin等(2014) [10] 研究了一步Newton-Raphson估计和局部轮廓似然估计，并给出了基于两种方法的估计量的分布；Alharbi与Patili (2018) [11] 提出对响应变量和残差的乘积进行平滑处理，并研究了基于此方法得到的方差函数估计量的渐近分布；Li和Lin (2020) [12] 在没有独立性假设的情形下，推导出了误差方差 ${\sigma }^2=\mathrm{var}\left(\epsilon \right)$ 的最佳半参数效率约束，并建立了基于残差的 ${\sigma }^2$ 有效估计量的渐近正态性。

本文对非参数回归模型中基于样条方法的均值函数估计量和方差函数估计量的分布问题进行研究，并通过数值模拟验证效果。

2. 估计量及主要结论

根据文献 [7]，将区间 $D=\left[0,1\right]$ 进行包括两端点的 $k+1$ 等距分割，结点序列为

$0=t_0<t_1<\cdots <t_{k+1}=1$，

构造相应的v次样条空间 $S_{k,v}$，其基函数记作 $B_s\left(x\right)\left(s=1,2,\cdots ,k+v\right)$。设 $K=k+v$，又令

$B_s\left(x\right)={\left(B_1\left(x\right),B_2\left(x\right),\cdots ,B_K\left(x\right)\right)}^{\text{T}}$，

则均值函数 $m\left(x\right)$ 的样条估计为

$\hat{m}\left(x\right)=B^{\text{T}}\left(x\right)\hat{\phi }={\displaystyle \underset{s\text{=}1}{\overset{K}{\sum }}{\hat{\phi }}_s}{B}_s\left(x\right)$， (2)

这里 $\hat{\phi }={\left({\hat{\phi }}_1,\cdots ,{\hat{\phi }}_K\right)}^{\text{T}}$ 为能使

$l_m\left(\phi \right)={\displaystyle \underset{t=1}{\overset{n}{\sum }}{\left\{y_t-{\displaystyle \underset{s=1}{\overset{K}{\sum }}{\phi }_s{B}_s\left(x_t\right)}\right\}}^2}$

最小化的参数向量， $\phi ={\left({\phi }_1,\cdots ,{\phi }_K\right)}^{\text{T}}$，它的最小二乘估计为 $\hat{\phi }={\left(W^{\text{T}}W\right)}^{-1}{W}^{\text{T}}y$，其中

$y={\left(y_1,y_2,\cdots ,y_n\right)}^{\text{T}}$，

$W=\left(\begin{array}{cccc}{B}_{k1}\left(x_1\right)& B_{k2}\left(x_1\right)& \cdots & B_{kK}\left(x_1\right)\\ B_{k1}\left(x_2\right)& B_{k2}\left(x_2\right)& \cdots & B_{kK}\left(x_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ B_{k1}\left(x_n\right)& B_{k2}\left(x_n\right)& \cdots & B_{kK}\left(x_n\right)\end{array}\right)$。

记残差 ${\hat{\varsigma }}_t=y_t-\hat{m}\left(x_t\right)\text{,}t=1,\cdots ,n$，令 $z_t={\hat{\varsigma }}_t^2$， $z={\left(z_1,\cdots ,z_n\right)}^{\text{T}}$，则方差函数 ${\sigma }^2\left(x\right)$ 的基于残差的样条估计为

${\hat{\sigma }}^2\left(x\right)=B^{\text{T}}\left(x\right)\hat{\theta }={\displaystyle \underset{s=1}{\overset{K}{\sum }}{\hat{\theta }}_s{B}_s\left(x\right)}$， (3)

其中 $\hat{\theta }={\left({\hat{\theta }}_1,\cdots ,{\hat{\theta }}_K\right)}^{\text{T}}$ 为使得

$l_{{\sigma }^2}\left(\theta \right)={{\displaystyle \underset{t=1}{\overset{n}{\sum }}\left\{z_t-{\displaystyle \underset{s=1}{\overset{K}{\sum }}{\theta }_s{B}_s\left(x_t\right)}\right\}}}^2$

最小化的参数向量， $\theta ={\left({\theta }_1,\cdots ,{\theta }_K\right)}^{\text{T}}$，它的最小二乘估计为 $\hat{\theta }={\left(W^{\text{T}}W\right)}^{-1}{W}^{\text{T}}z$。

2.1. 均值函数估计量的分布

本文假定误差序列 $\left\{{\epsilon }_t\right\}$ 来自标准正态总体，即 ${\epsilon }_t\sim N\left(0,1\right),t=1,2,\cdots ,n$。不妨记

$M={\left(m\left(x_1\right),m\left(x_2\right),\cdots ,m\left(x_n\right)\right)}^{\text{T}}$，

$\sum =\left(\begin{array}{cccc}\sigma \left(x_1\right)& 0& \cdots & 0\\ 0& \sigma \left(x_2\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \sigma \left(x_n\right)\end{array}\right)$，

$\epsilon ={\left({\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_n\right)}^{\text{T}}$，

则模型(1)可简记为

$y=M+\sum \epsilon$， (4)

其中

$\epsilon ~N\left(\left(\begin{array}{c}0\\ 0\\ ⋮\\ 0\end{array}\right),\left(\begin{array}{cccc}1& 0& \cdots & 0\\ 0& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\end{array}\right)\right)$。

因为 $M$ 与 $\sum$ 分别为常数向量和常数矩阵，又根据期望与方差的性质可知， $y~N\left(M,{\sum }^2\right)$。

定理1均值函数 $m\left(x\right)$ 的样条估计量 $\hat{m}\left(x\right)$ 服从均值为 $B^{\text{T}}\left(x\right)AM$，方差为 $B^{\text{T}}\left(x\right)A{\sum }^2{A}^{\text{T}}B\left(x\right)$ 的正态分布，即

$\hat{m}\left(x\right)~N\left(B^{\text{T}}\left(x\right)AM,B^{\text{T}}\left(x\right)A{\sum }^2{A}^{\text{T}}B\left(x\right)\right)$，

其中 $A={\left(W^{\text{T}}W\right)}^{-1}{W}^{\text{T}}$。

证明：根据(2)式知，要证 $\hat{m}\left(x\right)$ 的分布，只需证 $\hat{\phi }$ 的分布，又因为 $\hat{\phi }={\left(W^{\text{T}}W\right)}^{-1}{W}^{\text{T}}y$，那么有

$\hat{\phi }~N\left({\left(W^{\text{T}}W\right)}^{-1}{W}^{\text{T}}M,{\left(W^{\text{T}}W\right)}^{-1}{W}^{\text{T}}{\sum }^{2}W{\left(W^{\text{T}}W\right)}^{-1}\right)$，

故有

$\hat{m}\left(x\right)=B^{\text{T}}\left(x\right)\hat{\phi }~N\left(B^{\text{T}}\left(x\right)AM,B^{\text{T}}\left(x\right)A{\sum }^2{A}^{\text{T}}B\left(x\right)\right)$。

2.2. 方差函数估计量的分布

本文记 $X\leftrightarrow Y$ 表示随机变量序列X与Y渐近同分布。接下来讨论方差函数估计量 ${\hat{\sigma }}^2\left(x\right)$ 的渐近分布。

定理2方差函数 ${\sigma }^2\left(x\right)$ 基于残差的样条估计量 ${\hat{\sigma }}^2\left(x\right)$ 与 ${\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t{\epsilon }_t^2}$ 渐近同分布，即

${\hat{\sigma }}^2\left(x\right)\leftrightarrow {\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t{\epsilon }_t^2}$，

这里 $c=\left(c_1,c_2,\cdots ,c_n\right)=B^{\text{T}}\left(x\right)A{\sum }^2$。

证明：根据(3)式，要证方差函数估计量 ${\hat{\sigma }}^2\left(x\right)$ 的分布，只需证 $z={\left(z_1,z_2,\cdots ,z_n\right)}^{\text{T}}$ 的分布。

因为

${\hat{\varsigma }}_t=y_t-\hat{m}\left(x_t\right)=y_t-m\left(x_t\right)+m\left(x_t\right)-\hat{m}\left(x_t\right)=\sigma \left(x_t\right){\epsilon }_t-\left[\hat{m}\left(x_t\right)-m\left(x_t\right)\right]$，

则有

$z_t={\hat{\varsigma }}_t^2={\left\{\sigma \left(x_t\right){\epsilon }_t-\left[\hat{m}\left(x_t\right)-m\left(x_t\right)\right]\right\}}^2$，

即

$z_t={\left[\sigma \left(x_t\right){\epsilon }_t\right]}^2-2\sigma \left(x_t\right){\epsilon }_t\left[\hat{m}\left(x_t\right)-m\left(x_t\right)\right]+{\left[\hat{m}\left(x_t\right)-m\left(x_t\right)\right]}^2$。

由 $\hat{m}\left(x\right)$ 一致收敛到 $m\left(x\right)$ (见文献 [7] 中定理1)且 ${\epsilon }_t$ 是正态随机变量，可知

${\hat{\sigma }}^2\left(x\right)\leftrightarrow B^{\text{T}}\left(x\right)\cdot {\left(W^{\text{T}}W\right)}^{-1}{W}^{\text{T}}\cdot \left(\begin{array}{c}{\sigma }^2\left(x_1\right){\epsilon }_1^2\\ {\sigma }^2\left(x_2\right){\epsilon }_2^2\\ ⋮\\ {\sigma }^2\left(x_n\right){\epsilon }_n^2\end{array}\right)=B^{\text{T}}\left(x\right)A{\sum }^2\left(\begin{array}{c}{\epsilon }_1^2\\ {\epsilon }_2^2\\ ⋮\\ {\epsilon }_n^2\end{array}\right)$，

即

${\hat{\sigma }}^2\left(x\right)\leftrightarrow {\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t{\epsilon }_t^2}$。

因为 $\left\{{\epsilon }_t\right\}$ 为相互独立且服从标准正态分布的随机误差序列，故有

${\epsilon }_t^2~{\chi }^2\left(1\right),\forall t=1,2,\cdots ,n$，

因此，要求方差函数估计量 ${\hat{\sigma }}^2\left(x\right)$ 的渐近分布就是求服从卡方分布的独立随机变量的线性组合的分布。根据文献 [13] [14] 的结果，尝试用单个 ${\chi }^2$ 变量的线性函数近似n个相互独立的 ${\chi }^2$ 变量的线性组合。

首先，考虑用 $\tilde{Y}=a{\chi }^2\left(d\right)$ 作为 ${\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t{\epsilon }_t^2}$ 的近似分布。采用一、二阶矩拟合的原则选取 $a,d$，即由方程

$\{\begin{array}{l}E\left(\tilde{Y}\right)=E\left({\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t{\epsilon }_t^2}\right)\\ D\left(\tilde{Y}\right)=D\left({\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t{\epsilon }_t^2}\right)\end{array}$

确定。从而 $a,d$ 应满足方程

$\{\begin{array}{l}ad={\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t}\triangleq P\\ a^{2}d={\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t^2}\triangleq Q\end{array}$

解得

$a=Q/P,d=P^2/Q$。

考虑到 ${\chi }^2$ 变量的自由度为正整数，故将d修正为

$d^{*}=\left[P^2/Q+0.5\right]$，

这里 $\left[x\right]$ 表示不超过x的最大整数，若 $P^2/Q<0.5$，则取 $d^{*}=1$。

再用 $Y~a{\chi }^2\left(d^{*}\right)+e_1$ 来近似 ${\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t{\epsilon }_t^2}$ 的分布，采用上述方法得到

$e_1=P-\sqrt{Q{d}^{*}},a=\left(P-e_1\right)/d^{*}$，

于是，可用 $a{\chi }^2\left(d^{*}\right)+e_1$ 作为 ${\displaystyle \underset{t=1}{\overset{n}{\sum }}{c}_t{\epsilon }_t^2}$ 的近似。

3. 数值模拟

考虑模型(1)，其中

$\{\begin{array}{l}m\left(x\right)=50x^3{\left(1-x\right)}^3\\ \sigma \left(x\right)=0.2+0.4\sin \left(\pi x\right)，\\ {\epsilon }_t~N\left(0,1\right),t=1,2,\cdots ,n\end{array}$ (5)

这里 $x\in \left[0,1\right]$。

对模型(5)进行蒙特卡罗模拟，利用三次B-样条基函数估计未知均值函数 $m\left(x\right)$，基于残差估计方差函数 ${\sigma }^2\left(x\right)$，并在AIC准则下自动选取等距结点数。

3.1. 均值函数估计量的模拟

为验证理论分布效果，使用MATLAB软件进行模拟运算，选取显著性水平为 $\alpha =0.05$，具体步骤如下：

第一步，根据模型(5)，计算出 $B^{\text{T}}\left(x\right)AM$ 与 $B^{\text{T}}\left(x\right)A{\sum }^2{A}^{\text{T}}B\left(x\right)$ ；

第二步，生成N个服从 $N\left(B^{\text{T}}\left(x\right)AM,B^{\text{T}}\left(x\right)A{\sum }^2{A}^{\text{T}}B\left(x\right)\right)$ 的随机数；

第三步，基于样条方法，通过蒙特卡罗模拟产生N个均值函数的估计值 $\hat{m}\left(x\right)$ ；

第四步，在给定置信水平下，检验第二步与第三步产生的随机数与均值函数估计值是否来自于同一分布；

第五步，对上述四步进行多次重复模拟，分析所得结果。

图1绘制了 $x=0.02,0.2,0.4,0.6,0.8,1.0$ 处在 $N=1000$ 时的直方图和概率密度函数曲线图。

由图1可见，各点处拟合的均值函数的直方图与概率密度函数曲线呈倒U型，直观地可以认为各点处拟合的均值函数估计值来自于正态分布。进一步地，对各点处的均值函数估计量的分布与正态分布进行Two-sample t-test检验，并分别循环模拟 $N=10,50,100,500,1000$ 次，检验的P值如表1所示。

Figure 1. Histogram of the estimates at each point of the mean function and its probability density function curve

图1. 均值函数各点处估计值的直方图及其概率密度函数曲线

Table 1. The P-values of the Two-sample t-test-test for the mean function

表1. 均值函数Two-sample t-test检验的P值

由表1可知，当 $x=0.2,\text{0}\text{.4,}0.8$ 时，P值均大于0.05；当N较大时，各点处的P值均大于0.05，说明在给定的显著性水平0.05下，应该接受原假设，即认为检验数据服从正态分布。

3.2. 方差函数估计量的模拟

检验步骤如下：

第一步：依据模型(5)给定x值，计算a， $d^{*}$ 与 $e_1$ ；

第二步：随机生成N个服从 $a{\chi }^2\left(d^{*}\right)+e_1$ 的随机数；

第三步：基于残差样条方法，通过蒙特卡罗模拟生成N个方差函数 ${\hat{\sigma }}^2\left(x\right)$ 的估计值；

第四步，在给定置信水平下，检验第二步与第三步产生的随机数与方差函数估计值是否来自于同一分布；

第五步，对上述四步进行多次重复模拟，分析所得结果。

图2绘制了 $x=0.02,0.2,0.4,0.6,0.8,1.0$ 处在 $N=1000$ 时的直方图和概率密度函数曲线图。

Figure 2. Histogram of the estimates at each point and its probability density function curve

图2. 方差函数各点处估计值的直方图及其概率密度函数曲线

由图2可见，各点处拟合的方差函数估计量的直方图与概率密度函数曲线呈不对称的倒U型，且整体上右偏，直观地可以认为各点处拟合的方差函数估计量近似服从卡方分布。进一步地，对各点处的方差函数估计量的渐近分布与 $a{\chi }^2\left(d^{*}\right)+e_1$ 进行Two-sample t-test检验，并分别循环模拟 $N=10,50,100,500,1000$ 次，结果如表2所示。

由表2可知，各点处的P值均大于给定的显著性水平0.05，说明可以接受原假设，即认为检验数据近似服从 $a{\chi }^2\left(d^{*}\right)+e_1$。

Table 2. The P-values of the Two-sample t-test-test for the variance function

表2. 方差函数Two-sample t-test检验的P值

4. 结论

本文基于样条方法研究了固定设计下异方差非参数回归模型的均值函数与方差函数估计量的分布，均值函数的估计量服从正态分布，方差函数估计量的渐近分布可由单个 ${\chi }^2$ 变量的线性函数来近似。模拟结果显示：在给定显著性水平0.05下，分布拟合效果较优。

本文所研究的固定设计下异方差非参数回归模型的均值函数与方差函数估计量的近似分布为生物、医学、地质、经济等领域的研究带来了便利。

基金项目

国家自然科学基金项目(11601126)；河南省重点攻关项目(182102210286)。

NOTES

^*通讯作者。