空间变系数模型的局部非线性BGWR估计

doi:10.12677/AAM.2018.710144

期刊菜单

空间变系数模型的局部非线性BGWR估计
Local Nonlinear BGWR Estimation of Spatially Varying Coefficient Model

DOI: 10.12677/AAM.2018.710144, PDF, HTML, XML,
作者: 张晨晨, 张辉国, 胡锡健^*：新疆大学数学与系统科学学院，新疆乌鲁木齐
关键词: 空间变系数模型；局部非线性GWR；局部非线性BGWR；Gibbs抽样；Spatially Varying Coefficient Model； Local Nonlinear GWR； Local Nonlinear BGWR； Gibbs Sampling

摘要: 基于空间变系数模型的局部非线性GWR拟合方法的基础上，给出了空间变系数模型的局部非线性BGWR拟合方法。通过Gibbs抽样得到模型参数在各个地理位置处的估计值。同时，设计模拟实验，通过可视化形式，将该方法与局部非线性GWR估计方法作对比，显示出局部非线性BGWR方法的精确性。

Abstract: Based on the local nonlinear GWR fitting method of the spatial variable coefficient model, the local nonlinear BGWR fitting method of the spatial variable coefficient model is given. Estimates of model parameters at various geographic locations are obtained by Gibbs sampling. At the same time, the simulation experiment is designed. By comparing the method with the local nonlinear GWR estimation method, the accuracy of the local nonlinear BGWR method is shown.

文章引用：张晨晨, 张辉国, 胡锡健. 空间变系数模型的局部非线性BGWR估计[J]. 应用数学进展, 2018, 7(10): 1241-1246. https://doi.org/10.12677/AAM.2018.710144

1. 引言

空间变系数模型是线性模型适应空间数据的进一步推广，它将数据的地理位置引入回归系数中，以系数函数在空间位置处的估计值来分析回归关系的空间非平稳性特征，该模型满足了人们对处理空间问题的需要。

Brunsdon等 [1] 首次明确提出了空间变系数模型，同时给出了该模型的地理加权回归(geographical weighed regression，简称GWR)方法，即以观测点之间的距离函数为权值的非参数光滑估计方法。Wang Ning等 [2] [3] 将变系数模型的局部线性拟合方法推广到空间变系数模型的情况，给出了局部线性GWR方法。程鹏鹏 [4] 提出了局部非线性GWR拟合方法，在降低系数函数估计的偏和边界效应上较局部线性GWR方法更显著。LeSage [5] [6] 提出了贝叶斯地理加权回归(Bayesian geographical weighed regression，简称BGWR)模型，首先运用Gibbs抽样技术给出一组参数用来修正异常观测值，由贝叶斯公式得到待估参数的条件后验概率 [7] ，保证估计结果的稳健性，之后再根据GWR方法的原理进行估计。冯金杰 [8] 通过空间加权距离构造权重矩阵，基于多元线性回归模型的贝叶斯统计推断 [9] ，得到了该模型的局部线性BGWR估计方法。

本文在局部线性BGWR拟合方法的基础上，给出空间变系数模型的局部非线性BGWR估计方法，推导出各回归系数的满条件后验分布，采取Gibbs抽样方法得到模型参数的估计。通过模拟，将此方法的参数拟合结果可视化、计算偏差均值和标准差均值，并与局部非线性GWR估计方法的结果进行对比分析。

2. 空间变系数模型的局部非线性BGWR估计

本文的空间变系数自回归模型为：

$Y_{i} = \sum_{j = 1}^{p} β_{j} (u_{i}, v_{i}) X_{i j} + ε_{i}, i = 1, 2, \dots, n$ (1)

令模型中的 $X_{i 1} \equiv 1$ 是模型所包含得空间截距项， $Y = {(Y_{1}, Y_{2}, \dots, Y_{n})}^{T}$ ， $ε ~ N_{n} (0, σ^{2} I_{n})$ $ε = {(ε_{1}, ε_{2}, \dots, ε_{n})}^{T}$ ，由泰勒公式知在 $(u_{0}, v_{0})$ 的邻域内 $β_{j} (u, v)$ 可进行二阶泰勒展开，并用u和v的二次函数作逼近，即有

$\begin{array}{l} β_{j} (u, v) \approx β_{j} (u_{0}, v_{0}) + β_{j}^{(u)} (u_{0}, v_{0}) (u - u_{0}) + β_{j}^{(v)} (u_{0}, v_{0}) (v - v_{0}) + \frac{1}{2!} β_{j}^{(u u)} (u_{0}, v_{0}) {(u - u_{0})}^{2} \\ + \frac{1}{2!} β_{j}^{(v v)} (u_{0}, v_{0}) {(v - v_{0})}^{2} + \frac{1}{2!} β_{j}^{(u v)} (u_{0}, v_{0}) (u - u_{0}) (v - v_{0}) \end{array}$ (2)

其中 $β_{j}^{(u)} (u_{0}, v_{0})$ ， $β_{j}^{(v)} (u_{0}, v_{0})$ ， $β_{j}^{(u u)} (u_{0}, v_{0})$ ， $β_{j}^{(v v)} (u_{0}, v_{0})$ ， $β_{j}^{(u v)} (u_{0}, v_{0})$ 分别表示 $β_{j} (u, v)$ 关于u和v的一阶、二阶偏导数在 $(u_{0}, v_{0})$ 处的值。根据变系数模型的局部非线性拟合及GWR方法，选取未知参数 $β_{j} (u_{0}, v_{0})$ ， $β_{j}^{(v)} (u_{0}, v_{0})$ ， $β_{j}^{(u u)} (u_{0}, v_{0})$ ， $β_{j}^{(v v)} (u_{0}, v_{0})$ 和 $β_{j}^{(u v)} (u_{0}, v_{0})$ $(j = 1, 2, \dots, q)$ ，使得：

$\begin{array}{l} \sum_{i = 1}^{n} [Y_{i} - \sum_{j = 1}^{p} (β_{j} (u_{0}, v_{0}) + β_{j}^{(u)} (u_{0}, v_{0}) (u - u_{0}) + β_{j}^{(v)} (u_{0}, v_{0}) (v - v_{0}) + \frac{1}{2!} β_{j}^{(u u)} (u_{0}, v_{0}) {(u - u_{0})}^{2} \\ + \frac{1}{2!} β_{j}^{(v v)} (u_{0}, v_{0}) {(v - v_{0})}^{2} + β_{j}^{(u v)} (u_{0}, v_{0}) (u - u_{0}) (v - v_{0})) {X_{i j}]}^{2} K_{h} (d_{0 i}) \end{array}$ (3)

达到最小，其中 $(Y_{i}, X_{i 1}, X_{i 2}, \dots, X_{i p - 1}, X_{i p})$ 分别为因变量Y和自变量 $X_{1}, X_{2}, \dots, X_{p - 1}, X_{p}$ 在地理位置 $(u_{i}, v_{i})$ 处的观测值， $d_{0 i}$ 为 $(u_{0}, v_{0})$ 到 $(u_{i}, v_{i})$ 的距离，

$K_{h} (d_{0 i}) = K (\frac{d_{0 i}}{h}) = ω_{i} (u_{0}, v_{0}) (i = 1, 2, \dots, n)$

$W (u_{0}, v_{0}) = D i a g (ω_{1} (u_{0}, v_{0}), ω_{2} (u_{0}, v_{0}), \dots, ω_{n} (u_{0}, v_{0}))$

令

$\begin{matrix} P (u_{0}, v_{0}) = [β_{1} (u_{0}, v_{0}), \dots, β_{p} (u_{0}, v_{0}), β_{1}^{(u)} (u_{0}, v_{0}), \dots, β_{p}^{(u)} (u_{0}, v_{0}), \\ β_{1}^{(v)} (u_{0}, v_{0}), \dots, β_{p}^{(v)} (u_{0}, v_{0}), \dots, β_{1}^{(u u)} (u_{0}, v_{0}), \dots, β_{p}^{(u u)} (u_{0}, v_{0}), \\ β_{1}^{(v v)} (u_{0}, v_{0}), \dots, β_{p}^{(v v)} (u_{0}, v_{0}), β_{1}^{(u v)} (u_{0}, v_{0}), \dots, {β_{p}^{(u v)} (u_{0}, v_{0})]}^{T} \end{matrix}$

设计矩阵为 $X (u_{0}, v_{0}) = (A_{1}, A_{2}, A_{3}, A_{4})$ ，其中：

$A_{1} = (\begin{matrix} X_{11} & \dots & X_{1, p - 1} & X_{1 p} \\ X_{21} & \dots & X_{2, p - 1} & X_{2 p} \\ ⋮ & ⋮ & ⋮ \\ X_{n 1} & \dots & X_{n, p - 1} & X_{n p} \end{matrix})$

$A_{2} = (\begin{matrix} X_{11} (u_{1} - u_{0}) & \dots & X_{1, p - 1} (u_{1} - u_{0}) & X_{1 p} (u_{1} - u_{0}) \\ X_{21} (u_{2} - u_{0}) & \dots & X_{2, p - 1} (u_{2} - u_{0}) & X_{2 p} (u_{2} - u_{0}) \\ ⋮ & ⋮ & ⋮ \\ X_{n 1} (u_{n} - u_{0}) & \dots & X_{n, p - 1} (u_{n} - u_{0}) & X_{n p} (u_{n} - u_{0}) \end{matrix} \begin{matrix} X_{11} (v_{1} - v_{0}) & \dots & X_{1, p - 1} (v_{1} - v_{0}) & X_{1 p} (v_{1} - v_{0}) \\ X_{21} (v_{2} - v_{0}) & \dots & X_{2, p - 1} (v_{2} - v_{0}) & X_{2 p} (v_{2} - v_{0}) \\ ⋮ & ⋮ & ⋮ \\ X_{n 1} (v_{n} - v_{0}) & \dots & X_{n, p - 1} (v_{n} - v_{0}) & X_{n p} (v_{n} - v_{0}) \end{matrix})$

$A_{3} = (\begin{matrix} X_{11} {(u_{1} - u_{0})}^{2} & \dots & X_{1, p - 1} {(u_{1} - u_{0})}^{2} & X_{1 p} {(u_{1} - u_{0})}^{2} \\ X_{21} {(u_{2} - u_{0})}^{2} & \dots & X_{2, p - 1} {(u_{2} - u_{0})}^{2} & X_{2 p} {(u_{2} - u_{0})}^{2} \\ ⋮ & ⋮ & ⋮ \\ X_{n 1} {(u_{n} - u_{0})}^{2} & \dots & X_{n, p - 1} {(u_{n} - u_{0})}^{2} & X_{n p} {(u_{n} - u_{0})}^{2} \end{matrix} \begin{matrix} X_{11} {(v_{1} - v_{0})}^{2} & \dots & X_{1, p - 1} {(v_{1} - v_{0})}^{2} & X_{1 p} {(v_{1} - v_{0})}^{2} \\ X_{21} {(v_{2} - v_{0})}^{2} & \dots & X_{2, p - 1} {(v_{2} - v_{0})}^{2} & X_{2 p} {(v_{2} - v_{0})}^{2} \\ ⋮ & ⋮ & ⋮ \\ X_{n 1} {(v_{n} - v_{0})}^{2} & \dots & X_{n, p - 1} {(v_{n} - v_{0})}^{2} & X_{n p} {(v_{n} - v_{0})}^{2} \end{matrix})$

$A_{4} = (\begin{matrix} X_{11} (u_{1} - u_{0}) (v_{1} - v_{0}) & \dots & X_{1, p - 1} (u_{1} - u_{0}) (v_{1} - v_{0}) & X_{1 p} (u_{1} - u_{0}) (v_{1} - v_{0}) \\ X_{21} (u_{2} - u_{0}) (v_{2} - v_{0}) & \dots & X_{2, p - 1} (u_{2} - u_{0}) (v_{2} - v_{0}) & X_{2 p} (u_{2} - u_{0}) (v_{2} - v_{0}) \\ ⋮ & ⋮ & ⋮ \\ X_{n 1} (u_{n} - u_{0}) (v_{n} - v_{0}) & \dots & X_{n, p - 1} (u_{n} - u_{0}) (v_{n} - v_{0}) & X_{n p} (u_{n} - u_{0}) (v_{n} - v_{0}) \end{matrix})$

则局部非线性GWR方法的系数函数估计值为：

$\hat{B} (u_{0}, v_{0}) = {[X^{T} (u_{0}, v_{0}) W (u_{0}, v_{0}) X (u_{0}, v_{0})]}^{- 1} X^{T} (u_{0}, v_{0}) W (u_{0}, v_{0}) Y$ (4)

于是有，，因此似然函数为：

$\begin{matrix} L [B (u_{0}, v_{0}), σ] = {(\frac{1}{2 π σ^{2}})}^{\frac{n}{2}} \exp {- \frac{1}{2 σ^{2}} {[Y - X (u_{0}, v_{0}) B (u_{0}, v_{0})]}^{T} [Y - X (u_{0}, v_{0}) B (u_{0}, v_{0})]} \\ = {(\frac{1}{2 π σ^{2}})}^{\frac{n}{2}} \exp {- \frac{1}{2 σ^{2}} [{[Y - X (u_{0}, v_{0}) \hat{B} (u_{0}, v_{0})]}^{T} [Y - X (u_{0}, v_{0}) \hat{B} (u_{0}, v_{0})] \\ _{}^{} + {[B (u_{0}, v_{0}) - \hat{B} (u_{0}, v_{0})]}^{T} X^{T} (u_{0}, v_{0}) X (u_{0}, v_{0}) [B (u_{0}, v_{0}) - \hat{B} (u_{0}, v_{0})]]} \\ = {(\frac{1}{2 π σ^{2}})}^{\frac{n}{2}} \exp {- \frac{1}{2 σ^{2}} [Z_{n}^{2} (n - 6 p) + {[B (u_{0}, v_{0}) - \hat{B} (u_{0}, v_{0})]}^{T} \\ _{}^{} \cdot X^{T} (u_{0}, v_{0}) X (u_{0}, v_{0}) [B (u_{0}, v_{0}) - \hat{B} (u_{0}, v_{0})]]} \end{matrix}$

其中 $Z_{n}^{2} = {[Y - X (u_{0}, v_{0}) \hat{B} (u_{0}, v_{0})]}^{T} [Y - X (u_{0}, v_{0}) \hat{B} (u_{0}, v_{0})] / (n - 6 p)$ 。

若取 $π [B (u_{0}, v_{0}) | σ] \propto 1$ ， $π (σ) \propto \frac{1}{σ}$ ，那么由贝叶斯定理可知 $π (B (u_{0}, v_{0}), σ) \propto \frac{1}{σ}$ ， $σ > 0$ ，且参数 $β$ 的满条件后验分布为：

$\begin{matrix} π (β (u_{0}, v_{0}) | Y, X) \propto L (β (u_{0}, v_{0}), σ | Y, X) π (β (u_{0}, v_{0}), σ) \\ \propto M_{t_{3 p}} (v, \hat{β} (u_{0}, v_{0}), Q_{0}) \end{matrix}$ (6)

其中自由度为 $v = n - 6 p$ ，位置参数为 $\hat{β} (u_{0}, v_{0})$ ，精度矩阵为，

$Z_{n}^{2} = {[Y - X (u_{0}, v_{0}) \hat{β} (u_{0}, v_{0})]}^{T} [Y - X (u_{0}, v_{0}) \hat{β} (u_{0}, v_{0})] / (n - 6 p) .$

方差 $σ^{2}$ 的满条件后验分布为：

$π (σ^{2} | Y, X) \propto \frac{1}{σ^{n + 2}} \exp {- \frac{(n - 6 p) Z_{n}^{2}}{2 σ^{2}}} \propto I G a m m a (n + 1, v Z_{n}^{2} / 2)$ (7)

其中形状参数为 $n + 1$ ，尺度参数是 $v Z_{n}^{2} / 2$ 。

3. 模拟试验

本模拟实验，指定一个正方形区域，以此区域左下角为坐标原点建立直角坐标系，将横、纵坐标进行 $m - 1$ 等分处理，以坐标系中每个点作为地理位置 $(u_{i}, v_{i})$ ，则样本容量为 $n = m^{2}$ 。我们取 $m = 25$ ，有 $n = 625$ ，在如上建立的直角坐标系下，观测点地理位置坐标可表示为：

$(u_{i}, v_{i}) = (0.5 \mod (\frac{i - 1}{25}), 0.5 int (\frac{i - 1}{25})), i = 1, 2, \dots, n$

其中 $\mod ((i - 1) / m)$ 表示 $i - 1$ 除以m的余数， $int ((i - 1) / m)$ 表示 $i - 1$ 除以m的商的整数部分。

从而用以产生数据的空间变系数模型为：

$Y_{i} = β_{1} (u_{i}, v_{i}) + β_{2} (u_{i}, v_{i}) X_{i} + ε_{i}, i = 1, 2, \dots, n$

其中，自变量观测值 $X_{1}, X_{2}$ 来自均匀分布 $U (0, 2)$ 中的n个随机数，误差 $ε_{1}, ε_{2}, \dots, ε_{n}$ 来自正态分布 $N (0, 0.25)$ 中的n个随机数。选取系数函数 $β_{1} (u, v)$ 和 $β_{2} (u, v)$ ：

$β_{1} (u, v) = \frac{1}{6} (u + v), β_{2} (u, v) = \frac{1}{3} u$

在模拟实验中，选择高斯函数作为核函数，光滑参数h可使用CV准则来确定。根据已推出的参数的满条件后验分布，通过Gibbs抽样实现参数的 $β (u_{0}, v_{0})$ 和 $σ^{2}$ 的估计，让 $(u_{0}, v_{0})$ 遍历所有位置，可得到所有点的系数估计值。

绘制出 $β_{1}$ ， $β_{2}$ 的初值曲面图(如图1)和各回归系数函数在不同参数估计方法下的估计曲面图(如图2a，2b)。相比之下，局部非线性BGWR估计方法得到的系数曲面更接近原曲面的基本特征，而局部非线性GWR估计得到的曲面边界较为粗糙。因此，在对模型系数函数估计时，局部非线性BGWR估计方法更具精确性。

分别计算两种估计方法的绝对偏差MBias和均方误差MSE的值，结果如表1。不难发现，局部非线性BGWR方法估计相比局部非线性GWR方法的 $M B i a s ({\hat{β}}_{j})$ 和 $M S E ({\hat{β}}_{j})$ 的值要小一些，说明局部非线性BGWR方法估计系数函数估计的精确度相比较高。

Figure 1. Surface map of the initial values of each coefficient function

图1. 系数函数初值曲面图

(a) (b)

Figure 2. (a) Estimation of each coefficient function of the local nonlinear GWR fitting method, (b) Estimation of each coefficient function of the local nonlinear BGWR fitting method

图2. (a) 局部非线性GWR拟合方法的模型系数函数估计值，(b) 局部非线性BGWR拟合方法的模型系数函数估计值

Table 1. Absolute deviation and mean square error of the estimated

表1. 估计值MBias和MSE结果

4. 结论

本文基于对空间变系数模型的参数拟合方法研究的思想，结合局部非线性GWR方法和贝叶斯方法，提出局部非线性BGWR拟合方法，推导出各参数的满条件后验分布，采用Gibbs抽样技术得到模型的系数函数估计值。同时，本文设计模拟实验，通过可视化形式，展现出局部非线性BGWR和局部非线性GWR方法的参数拟合效果，结果显示，采用局部非线性BGWR方法估计的各估计值的绝对偏差、均方误差值和系数函数拟合图均优于局部非线性GWR参数拟合方法。

NOTES

^*第一作者。

^#通讯作者。

参考文献

[1]	Brunsdon, C., Fotheringham, A.S. and Charlton, M. (1996) Geographically Weighted Regression: A Method for Exploring Spatial Non-Stationarity. Geographically Analysis, 28, 281-298. https://doi.org/10.1111/j.1538-4632.1996.tb00936.x
[2]	梅长林, 王宁. 近代回归分析方法[M]. 北京: 科学出版社, 2012.
[3]	Wang, N., Mei, C.L. and Yan, X.D. (2008) Local Linear Estimation of Spatially Varying Coefficient Models: An Improvement on the Geographically Weighted Regression Technique. Environment and Planning A, 40, 986-1005. https://doi.org/10.1068/a3941
[4]	程鹏鹏. 空间变系数模型的局部线性估计[J]. 哈尔滨师范大学自然科学学报, 2014, 30(3): 40-42.
[5]	LeSage, J.P. (2001) A Family of Geographically Weighted Regression Models. http://www.spatial-econometrics.com/html/bgwr.pdf
[6]	LeSage, J.P. (1997) Bayesian Estimation of Spatial Autoregressive Models. International Science Review, 20, 113-129. https://doi.org/10.1177/016001769702000107
[7]	詹姆斯·勒沙杰, 凯利·佩斯. 空间计量经济学导论[M]. 北京: 北京大学出版社, 2014.
[8]	冯金杰. 空间变系数模型的局部线性BGWR估计及其应用[D]: [硕士学位论文]. 乌鲁木齐: 新疆大学, 2016.
[9]	朱慧明, 韩玉启. 贝叶斯多元统计推断理论[M]. 北京: 科学出版社, 2006.

为你推荐

友情链接