浙江农林大学统计系，浙江 杭州

收稿日期：2020年11月26日；录用日期：2020年12月11日；发布日期：2020年12月18日

摘要

基于方差建模研究了变系数异方差模型的贝叶斯估计和异常点识别，其中非参数部分采用B样条逼近。主要通过应用Gibbs抽样和Metropolis-Hastings算法相结合的混合算法获得模型的贝叶斯估计和通过K-L距离贝叶斯诊断统计量来识别数据异常点。模拟研究显示所提出的贝叶斯分析方法是可行有效的。

关键词

异方差模型，Metropolis-Hastings算法，贝叶斯估计，K-L距离，B样条

Bayesian Analysis of Varying Coefficient Heteroscedastic Models

Fangzhong Xu, Dengke Xu

Department of Statistics, Zhejiang Agriculture and Forestry University, Hangzhou Zhejiang

Received: Nov. 26^th, 2020; accepted: Dec. 11^th, 2020; published: Dec. 18^th, 2020

ABSTRACT

Based on variance modeling, Bayesian estimation and outlier identification of varying coefficient heteroscedastic models are studied, where the nonparametric part is approximated by B-spline. By combining the Gibbs sampler and Metropolis-Hastings algorithm, Bayesian estimation and Bayesian diagnosis statistics based on the K-L distance are obtained to identify outliers. Simulation studies show that the proposed Bayesian methods are feasible and effective.

Keywords:Heteroscedastic Models, Metropolis-Hastings Algorithm, Bayesian Estimation, K-L Distance, B-Spline

Copyright © 2020 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

异方差数据常常出现在经济学、生物学、环境科学等领域，是很多统计学家研究的热点方向之一。其中处理此类数据最常用的方法是方差建模法，即不仅对均值建立回归模型，同时也对方差建立回归模型进行分析，有些文献称之为联合均值与方差模型。这个模型主要体现了对方差的重视，它能更好地解释数据变化的原因和规律。特别最近这些年已经有很多学者对基于方差建模的异方差模型研究了模型的参数估计、变量选择以及异方差检验等统计推断。例如，吴刘仓等 [1] 对联合均值与方差模型提出一种同时对均值模型和方差模型的变量选择方法；李双双等 [2] 研究了联合均值与方差混合专家回归模型的参数估计问题；赵远英等 [3] 对响应变量带有不可忽略缺失数据的联合均值与方差模型的贝叶斯估计问题进行了研究；戴琳等 [4] 基于联合均值与方差模型研究了模型的参数估计与基于数据删除模型考虑了统计诊断问题。其它类似的相关研究还可以具体参见文献 [5] [6] [7] [8]。但是这些文献大多数都是基于参数异方差模型展开统计分析的，很少有文献和学者研究基于方差建模的变系数异方差模型的贝叶斯参数估计、异常点识别等统计推断问题。

因此本文针对方差建模的变系数异方差模型，应用Gibbs抽样和Metropolis-Hastings算法相结合的混合算法研究模型的贝叶斯估计和基于K-L距离研究贝叶斯数据删除影响诊断方法，以识别模型的异常点。

2. 模型与符号

2.1. 基于方差建模的变系数异方差模型

针对异方差数据和基于方差建模的思想，提出了如下变系数异方差模型：

$\{\begin{array}{l}{y}_i=x_i^{\text{T}}\beta +z_i^{\text{T}}\alpha \left(u_i\right)+{\epsilon }_i,{\epsilon }_i~N\left(0,{\sigma }_i^2\right),\\ {\sigma }_i^2=g\left(h_i^{\text{T}}\gamma \right),\\ i=1,2,\cdots ,n.\end{array}$ (1)

其中 $y_i$ 为响应变量且令 $Y={\left(y_1,y_2,\cdots ,y_n\right)}^{\text{T}}$ ； $x_i={\left(x_{i1},\cdots ,x_{ip}\right)}^{\text{T}}$，$z_i={\left(z_{i1},\cdots ,z_{iq}\right)}^{\text{T}}$ 是与 $y_i$ 的均值相关的解释变量的观测值，且 $\beta ={\left({\beta }_1,\cdots ,{\beta }_p\right)}^{\text{T}}$ 是均值模型中一个 $p\times 1$ 维未知回归参数向量， $\alpha \left(u_i\right)={\left({\alpha }_1\left(u_i\right),{\alpha }_2\left(u_i\right),\cdots ,{\alpha }_d\left(u_i\right)\right)}^{\text{T}}$ 是一个 $d\times 1$ 维的未知函数参数； $h_i={\left(h_{i1},\cdots ,h_{iq}\right)}^{\text{T}}$ 是与 $y_i$ 的方差相关的解释变量的观测值，且 $\gamma ={\left({\gamma }_1,\cdots ,{\gamma }_q\right)}^{\text{T}}$ 是方差模型中一个 $q\times 1$ 维的未知回归参数向量。另外， $g(\cdot )$ 是一个已知函数，为了模型的可识别性，一般假设 $g(\cdot )$ 是一个单调函数且 $g(\cdot )>0$。在本文中 $g\left(x\right)=\exp \left(x\right)$。

由模型(1)，可以获得如下似然函数

$L\left(\beta ,\gamma ,\alpha (\cdot )|Y,X,Z,H,U\right)={\left(2\pi \right)}^{-\frac{n}{2}}{\displaystyle \underset{i=1}{\overset{n}{\prod }}{\left({\sigma }_i^2\right)}^{-\frac{1}{2}}}\exp \left(-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{{\left(y_i-x_i^{\text{T}}\beta -z_i^{\text{T}}\alpha \left(u_i\right)\right)}^2}{{\sigma }_i^2}}\right)$ (2)

其中 $X={\left(x_1^{\text{T}},\cdots ,x_n^{\text{T}}\right)}^{\text{T}},Z={\left(z_1^{\text{T}},\cdots ,z_n^{\text{T}}\right)}^{\text{T}},H={\left(h_1^{\text{T}},\cdots ,h_n^{\text{T}}\right)}^{\text{T}},U={\left(u_1,\cdots ,u_n\right)}^{\text{T}}$。因为 $\alpha (\cdot )$ 是非参数，(2)式还不能直接进行优化。因此，首先用B样条来逼近非参数函数 $\alpha (\cdot )$。具体如下，令 $0=s_0<s_1<\cdots <s_{k_n}<s_{k_n+1}=1$ 是[0,1]区间上的一个剖分。用 $s_i$ 作为内节点，那么就有阶为M和维数为 $L=k_n+M$ 的正则化B样条基函数，这也形成了线性样条空间的一个基。节点选择一般是样条光滑估计中的一个重要方面。类似于文献 [9]，内节点的数目选取为 $n^{\frac{1}{5}}$ 的整数部分。这样，由 ${\pi }^{\text{T}}\left(u\right){\lambda }_k$ 逼近 ${\alpha }_k\left(u\right)$，其中 $\pi \left(u\right)={\left({\pi }_1\left(u\right),\cdots ,{\pi }_L\left(u\right)\right)}^{\text{T}}$ 是基函数向量和 ${\lambda }_k\in R^L,k=1,2,\cdots ,d$。利用这些符号，(1)式中的均值模型可以写成以下形式：

${\mu }_i=x_i^{\text{T}}\beta +B_i^{\text{T}}\lambda$

其中 $B_i=I_d\otimes \pi \left(u_i\right)\cdot z_i,\lambda ={\left({\lambda }_1^{\text{T}},\cdots ,{\lambda }_d^{\text{T}}\right)}^{\text{T}},K=dL$。

这样，似然函数(2)可以被重新写成以下形式：

$L\left(\theta |Y,X,Z,H,U\right)={\left(2\pi \right)}^{-\frac{n}{2}}{\displaystyle \underset{i=1}{\overset{n}{\prod }}{\left({\sigma }_i^2\right)}^{-\frac{1}{2}}}\exp \left(\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{{\left(y_i-x_i^{\text{T}}\beta -B_i^{\text{T}}\lambda \right)}^2}{{\sigma }_i^2}}\right)$，(3)

其中 $\theta ={\left({\beta }^{\text{T}},{\gamma }^{\text{T}},{\lambda }^{\text{T}}\right)}^{\text{T}}$。

2.2. K-L距离

K-L距离也称作K-L信息，具有距离和信息的某些性质，在统计上以反映两个模型或分布的差异而著称。根据韦博成 [10] 对K-L距离的介绍，密度函数为 $f\left(x\right)$ 与 $g\left(x\right)$ 的两个分布的K-L距离 $K\left(f,g\right)$ 被定义为： $K\left(f,g\right)=E_f\left\{\log \frac{f\left(x\right)}{g\left(x\right)}\right\}$，其中 $E_f$ 表示对 $f\left(x\right)$ 求期望。

3. 贝叶斯分析

3.1. 先验分布

为了应用贝叶斯方法来估计模型(1)中的未知参数，需要给出未知参数的先验分布信息。为了简便，假设 $\beta ,\gamma$ 和 $\lambda$ 相互独立且具有正态先验分布，分别为 $\beta ~N\left({\beta }_0,{\sum }_{\beta }\right)$，$\gamma ~N\left({\gamma }_0,{\sum }_{\gamma }\right)$，$\lambda ~N\left({\lambda }_0,{\tau }^2{I}_K\right)$，${\tau }^2~IG\left(a_{\tau },b_{\tau }\right)$，其中“IG”表示逆gamma分布，并且假设超参数 ${\beta }_0,{\gamma }_0,{\lambda }_0,{\sum }_{\beta },{\sum }_{\gamma }$ 和 $a_{\tau },b_{\tau }$ 是已知的。

3.2. Gibbs抽样和条件分布

基于式子(3)，按照以下过程用Gibbs抽样从后验分布 $p\left(\theta |Y,X,Z,H,U\right)$ 中进行抽样，其中 $\theta ={\left({\beta }^{\text{T}},{\gamma }^{\text{T}},{\lambda }^{\text{T}}\right)}^{\text{T}}$。

步骤 1. 令参数的初值 ${\theta }^{\left(0\right)}={\left({\beta }^{{\left(0\right)}^{\text{T}}},{\gamma }^{{\left(0\right)}^{\text{T}}},{\lambda }^{{\left(0\right)}^{\text{T}}}\right)}^{\text{T}}$。

步骤 2. 基于 ${\theta }^{\left(l\right)}={\left({\beta }^{{\left(l\right)}^{\text{T}}},{\gamma }^{{\left(l\right)}^{\text{T}}},{\lambda }^{{\left(l\right)}^{\text{T}}}\right)}^{\text{T}}$，计算 ${\sum }^{\left(l\right)}=diag\left\{{\sigma }_1^{2\left(l\right)},\cdots ,{\sigma }_n^{2\left(l\right)}\right\}$。

步骤 3. 基于 ${\theta }^{\left(l\right)}={\left({\beta }^{{\left(l\right)}^{\text{T}}},{\gamma }^{{\left(l\right)}^{\text{T}}},{\lambda }^{{\left(l\right)}^{\text{T}}}\right)}^{\text{T}}$ 按照以下抽取 ${\theta }^{\left(l+1\right)}={\left({\beta }^{{\left(l+1\right)}^{\text{T}}},{\gamma }^{{\left(l+1\right)}^{\text{T}}},{\lambda }^{{\left(l+1\right)}^{\text{T}}}\right)}^{\text{T}}$ ；

· 抽样 ${\tau }^{2\left(l+1\right)}$ ：

$p\left({\tau }^2|\lambda \right)\propto {\left({\tau }^2\right)}^{-\frac{K}{2}-a_{\tau }-1}\exp \left\{\frac{-\frac{1}{2}{\left(\lambda -{\lambda }_0\right)}^{\text{T}}\left(\lambda -{\lambda }_0\right)-b_{\tau }}{{\tau }^2}\right\}$ (4)

· 抽样 ${\lambda }^{\left(l+1\right)}$ ：

$p\left(\lambda |Y,X,Z,H,U,\beta ,\gamma \right)\propto \exp \left\{-\frac{1}{2}{\left(\lambda -{\lambda }^{*}\right)}^{\text{T}}{\Sigma }_{\lambda }^{\ast -1}\left(\lambda -{\lambda }^{*}\right)\right\}$， (5)

其中 ${\lambda }^{*}={\Sigma }_{\lambda }^{\ast }\left({\Omega }^{\text{T}}{\sum }^{-1}\left(Y-X\beta \right)+{\tau }^{-2}{I}_K{\lambda }_0\right)$，${\Sigma }_{\lambda }^{\ast }={\left({\tau }^{-2}{I}_K+{\Omega }^{\text{T}}{\Sigma }^{-1}\Omega \right)}^{-1}$，$\Omega ={\left(B_1^{\text{T}},\cdots ,B_n^{\text{T}}\right)}^{\text{T}}$。

· 抽样 ${\beta }^{\left(l+1\right)}$ ：

$p\left(\beta |Y,X,Z,H,U,\gamma ,\lambda \right)\propto \exp \left\{-\frac{1}{2}{\left(\beta -b^{*}\right)}^{\text{T}}{B}^{*-1}\left(\beta -b^{*}\right)\right\}$，(6)

其中 $b^{*}=B^{*}\left(X^T{\sum }^{-1}\left(Y-\Omega \lambda \right)+{\sum }_{\beta }^{-1}{\beta }_0\right)$，$B^{*}={\left({\sum }_{\beta }^{-1}+X^{\text{T}}{\Sigma }^{-1}X\right)}^{-1}$。

· 抽样 ${\gamma }^{\left(l+1\right)}$ ：

$p\left(\gamma |Y,X,Z,H,U,\beta ,\lambda \right)\propto \exp \left\{-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{h}_i^{\text{T}}\gamma -\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{{\left(y_i-x_i^{\text{T}}\beta -B_i^{\text{T}}\lambda \right)}^2}{\exp \left(h_i^{\text{T}}\gamma \right)}}-\frac{1}{2}{\left(\gamma -{\gamma }_0\right)}^{\text{T}}{\sum }_{\gamma }^{-1}\left(\gamma -{\gamma }_0\right)}\right\}$ (7)

步骤4。重复步骤2和3。

这样就通过以上算法产生了样本序列 $\left({\beta }^{\left(l\right)},{\gamma }^{\left(l\right)},{\lambda }^{\left(l\right)},{\tau }^2{}^{\left(l\right)}\right),l=1,2,\cdots$。从(4)~(7)式中很容易发现，条件分布 $p\left({\tau }^2|\lambda \right),p\left(\beta |Y,X,Z,H,U,\gamma ,\lambda \right),p\left(\lambda |Y,X,Z,H,U,\beta ,\gamma \right)$ 是熟悉的正态分布和逆gamma分布。从这两个分布抽取随机数是比较容易的。但是条件分布 $p\left(\gamma |Y,X,Z,H,U,\beta ,\lambda \right)$ 是一不规则且相当复杂的分布，如何从这个分布中抽取随机数有点困难。在这里主要应用MH算法抽取随机数。选择正态分布 $N\left({\gamma }^{\left(l\right)},{\sigma }_{\gamma }^2{\Omega }_{\gamma }^{-1}\right)$ 作为建议分布，其中通过选择 ${\sigma }_{\gamma }^2$，来使得接受概率在0.25与0.45之间，且取 ${\Omega }_{\gamma }={\sum }_{\gamma }^{-1}+\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{{\left(y_i-x_i^{\text{T}}\beta -B_i^{\text{T}}\lambda \right)}^2}{\exp \left(h_i^{\text{T}}\gamma \right)}{h}_i{h}_i^{\text{T}}}$。

3.3. 贝叶斯估计

利用以上提出的计算过程来产生观测值来获得参数 $\beta ,\gamma$ 和 $\lambda$ 的贝叶斯估计。令 $\left\{{\theta }^{\left(j\right)}=\left({\beta }^{\left(j\right)},{\gamma }^{\left(j\right)},{\lambda }^{\left(j\right)}\right):j=1,2,\cdots ,J\right\}$ 是通过上述混合算法从联合条件分布 $p\left(\beta ,\gamma ,\lambda |Y,X,Z,H,U\right)$ 中产生的观测值，那么 $\beta ,\gamma$ 和 $\lambda$ 的贝叶斯估计为：

$\hat{\beta }=\frac{1}{J}{\displaystyle \underset{j=1}{\overset{J}{\sum }}{\beta }^{\left(j\right)}},\hat{\gamma }=\frac{1}{J}{\displaystyle \underset{j=1}{\overset{J}{\sum }}{\gamma }^{\left(j\right)}},\hat{\lambda }=\frac{1}{J}{\displaystyle \underset{j=1}{\overset{J}{\sum }}{\lambda }^{\left(j\right)}}$。

根据文献 [11] 有，当J趋于无穷时， $\hat{\theta }=\left(\hat{\beta },\hat{\gamma },\hat{\lambda }\right)$ 是对应后验均值向量的相合估计。类似地，后验协方差矩阵 $Var\left(\theta |Y,X,Z,H,U\right)$ 的相合估计可以通过观测 $\left\{{\theta }^{\left(j\right)}:j=1,2,\cdots ,J\right\}$ 的样本协方差矩阵来获得，即

$Var\left(\theta |Y,X,Z,H,U\right)={\left(J-1\right)}^{-1}{\displaystyle \underset{j=1}{\overset{J}{\sum }}\left({\theta }^{\left(j\right)}-\hat{\theta }\right)}{\left({\theta }^{\left(j\right)}-\hat{\theta }\right)}^{\text{T}}$。

这样，后验标准误就可以通过该矩阵的对角元素来获得。

3.4. 贝叶斯诊断

在贝叶斯统计诊断分析中已存在许多诊断统计量用来评价个体观测对参数后验分布的影响，在这主要基于K-L距离研究贝叶斯数据删除影响的统计诊断方法。对任意的 $i=1,\cdots ,n$，记 $\left\{y_i,x_i,z_i,h_i,u_i\right\}$ 是第i个个体观测数据点， $D=\left\{Y,X,Z,H,U\right\}$ 为完全数据集， $D_{-i}$ 为完全数据集D删除第i个个体观测数据点得到的数据集， $L\left(\theta |D\right)$ 与 $L\left(\theta |D_{-i}\right)$ 分别表示基于数据D与 $D_{-i}$ 的似然函数，则 $\theta$ 于数据D与 $D_{-i}$ 的后验分布分别为 $p\left(\theta |D\right)\propto L\left(\theta |D\right)p\left(\theta \right)$，$p\left(\theta |D_{-i}\right)\propto L\left(\theta |D_{-i}\right)p\left(\theta \right)$。根据Cho等 [7] 的讨论，定义K-L距离为：

$K\left(P,P_{-i}\right)={\displaystyle \int p\left(\theta |D\right)\log \left\{\frac{p\left(\theta |D\right)}{p\left(\theta |D_{-i}\right)}\right\}\text{d}\theta }$， (8)

其中P与 $P_{-i}$ 分别表示 $\theta$ 基于数据D与 $D_{-i}$ 的后验分布，注意到K-L距离 $K\left(P,P_{-i}\right)$ 是完全数据集D删除第i个数据点前后对参数 $\theta$ 后验分布影响的一种很好的度量。经过简单的计算(8)式变为：

$K\left(P,P_{-i}\right)=\log E_{\theta }\left[\frac{L\left(\theta |D_{-i}\right)}{L\left(\theta |D\right)}|D\right]+E_{\theta }\left[\log \frac{L\left(\theta |D\right)}{L\left(\theta |D_{-i}\right)}|D\right]$， (9)

其中 $E_{\theta }\left[.|D\right]$ 表示 $\theta$ 基于数据D的后验期望。由Gibbs抽样算法抽取的随机观测序 $\left\{{\theta }^{\left(j\right)}:j=1,2,\cdots ,J\right\}$，可以得到K-L距离 $K\left(P,P_{-i}\right)$ 的估计为：

$K\left(P,P_{-i}\right)=\log \left[\frac{1}{J}{\displaystyle \underset{j=1}{\overset{J}{\sum }}\frac{L\left({\theta }^{\left(j\right)}|D_{-i}\right)}{L\left({\theta }^{\left(j\right)}|D\right)}}\right]+\frac{1}{J}{\displaystyle \underset{j=1}{\overset{J}{\sum }}\log \left[\frac{L\left({\theta }^{\left(j\right)}|D\right)}{L\left({\theta }^{\left(j\right)}|D_{-i}\right)}\right]}$。 (10)

对任意的 $i=1,\cdots ,n$，当 $K\left(P,P_{-i}\right)$ 很大时，可以诊断第i个个体观测数据点为异常点。

4. 模拟研究

在这通过模拟研究来说明本文所提出的贝叶斯分析方法的有效性。选择均值模型为 ${\mu }_i=x_i^{\text{T}}\beta +z_i^{\text{T}}\alpha \left(u_i\right)$ 和方差模型为 ${\sigma }_i^2=\exp \left(h_i^{\text{T}}\gamma \right)$，其中均值参数和方差参数的真实值分别为 $\beta ={\left(1,-0.5,0.5\right)}^{\text{T}}$ 和 $\gamma ={\left(1,-0.5,0.5\right)}^{\text{T}}$，另外 $d=2$，${\alpha }_1\left(u_i\right)=0.5\sin \left(2\pi u_i\right)$，${\alpha }_2\left(u_i\right)=8u_i\left(1-u_i^2\right)$ ； $x_i,z_i,h_i$ 分别是 $3\times 1,2\times 1$ 和 $3\times 1$ 的协变量向量，其中的元素产生于均值为零，协方差矩阵为 $\sum =\left({\sigma }_{ij}={0.5}^{\left|i-j\right|}\right)$ 的正态分布。这样响应变量 $Y_i$ 就可以从多元正态分布 $N\left({\mu }_i,{\sum }_i\right)\left(i=1,\cdots ,n\right)$ 中产生。

为了调查贝叶斯分析方法对先验分布的敏感程度，考虑以下三种有关未知参数 $\beta ,\gamma ,\lambda$ 的先验分布中超参数值的设置的情形：

TypeI：

${\beta }_0={\left(1,-0.5,0.5\right)}^{\text{T}},{\sum }_{\beta }=0.25\times {{\rm I}}_{\text{3}},{\gamma }_0={\left(1,-0.5,0.5\right)}^{\text{T}},{\sum }_{\gamma }=0.25\times {{\rm I}}_3,{\lambda }_0={\left(0,\cdots ,0\right)}^{\text{T}},a_{\tau }=b_{\tau }=1$。这种设置表示具有很好的先验信息。

TypeII：

${\beta }_0=3\times {\left(1,-0.5,0.5\right)}^{\text{T}},{\sum }_{\beta }={{\rm I}}_{\text{3}},{\gamma }_0=3\times {\left(1,-0.5,0.5\right)}^{\text{T}},{\sum }_{\gamma }={{\rm I}}_3,{\lambda }_0={\left(0,\cdots ,0\right)}^{\text{T}},a_{\tau }=b_{\tau }=1$。这种设置表示具有较差的先验信息。

TypeIII：

${\beta }_0={\left(0,0,0\right)}^{\text{T}},{\sum }_{\beta }=10\times {{\rm I}}_3,{\gamma }_0={\left(0,0,0\right)}^{\text{T}},{\sum }_{\gamma }=10\times {{\rm I}}_3,{\lambda }_0={\left(0,\cdots ,0\right)}^{\text{T}},a_{\tau }=b_{\tau }=1$。这些超参数值的设置代表的是没有先验信息的情况。

在上面的各种情形下，应用联合Gibbs抽样和Metropolis-Hastings算法的混合算法来计算未知参数和光滑函数的贝叶斯估计。在模拟中分别令样本量n = 80和n = 150。对于每一种情形，重复计算100次。对于每次重复产生的每一次数据集，MCMC算法的收敛性可以通过EPSR值来检验，并且在每次运行中观测得到在3000次迭代以后EPSR值都小于1.2。因此在每次重复计算中丢掉前3000次迭代以后再收集J = 5000个样本来产生贝叶斯估计。参数贝叶斯估计的模拟结果概括在表1中。为了调查估计函数 ${\alpha }_1\left(u\right)$ 和 ${\alpha }_2\left(u\right)$ 的精确度，画出了在不同样本量和三种先验分布情形下非参数函数的平均估计曲线和真实曲线，并且展示在图1~6中。

在表1中，“Bias”表示基于100次重复计算未知参数的贝叶斯估计与真值的偏差的绝对值，“SD”表示前面给出的后验标准误的平均估计，和“RMS”表示的是基于100次重复计算的贝叶斯估计的均方误差的算术平方根。从表1中可以获得i) 在估计的偏差、RMS和SD值方面，不管何种先验信息贝叶斯估计都相当精确；ii) 当样本量逐渐变大时，估计也变得越来越好。从图1~6中展示了不管何种先验信息，估计出来的非参数函数的曲线与相应的真实函数的曲线逼近得都比较好。总之，从以上结果可以看出本文所提出的贝叶斯估计方法能很好地恢复变系数异方差模型中的真实信息。

表1. 不同的样本量和先验分布下模型参数的贝叶斯估计结果

图1. 当n = 80和Type I先验信息下非参数部分 ${\alpha }_1\left(u\right)$ 和 ${\alpha }_2\left(u\right)$ 的平均估计曲线

图2. 当n = 80和Type II先验信息下非参数部分 ${\alpha }_1\left(u\right)$ 和 ${\alpha }_2\left(u\right)$ 的平均估计曲线

图3. 当n = 80和Type III先验信息下非参数部分 ${\alpha }_1\left(u\right)$ 和 ${\alpha }_2\left(u\right)$ 的平均估计曲线

图4. 当n = 150和Type I先验信息下非参数部分 ${\alpha }_1\left(u\right)$ 和 ${\alpha }_2\left(u\right)$ 的平均估计曲线

图5. 当n = 150和Type II先验信息下非参数部分 ${\alpha }_1\left(u\right)$ 和 ${\alpha }_2\left(u\right)$ 的平均估计曲线

图6. 当n = 150和Type III先验信息下非参数部分 ${\alpha }_1\left(u\right)$ 和 ${\alpha }_2\left(u\right)$ 的平均估计曲线

最后为了检测根据K-L距离来识别异常点的效果，在第10个个体观测数据点的响应变量加10构成人工数据集D。然后对人工数据集D应用本文介绍的贝叶斯影响诊断方法来检测影响观测。其中MCMC算法的收敛性可以通过EPSR值来检验，并且发现在3000次迭代以后EPSR值都小于1.2。因此在计算中丢掉前3000次迭代以后再收集J = 5000个随机样本来通过(10)式计算K-L距离 $K\left(P,P_{-i}\right)$。图7和图8报道了相应的诊断结果。正如预期的一样，通过图7和图8很容易就发现，第10个个体观测数据点被诊断为异常点，且诊断方法对先验分布超参数取值的选取不是很敏感。

图7. 当n = 80时，基于不同的先验信息下贝叶斯数据删除影响诊断的数值结果

图8. 当n = 150时，基于不同的先验信息下贝叶斯数据删除影响诊断的数值结果

5. 结论

本文针对变系数异方差模型，基于Gibbs抽样和MH算法相结合的混合算法，以及根据K-L距离研究模型的贝叶斯估计和贝叶斯统计诊断方法。模拟研究显示了模型与贝叶斯方法的可行性和有效性。

基金项目

浙江省高校重大人文社科攻关计划项目资助(2018QN037)。

文章引用

许芳忠,徐登可. 变系数异方差模型的贝叶斯分析
Bayesian Analysis of Varying Coefficient Heteroscedastic Models[J]. 应用数学进展, 2020, 09(12): 2166-2175. https://doi.org/10.12677/AAM.2020.912252

参考文献