DOI:10.12677/SA.2017.65064,PDF,HTML,XML,被引量下载: 1,677浏览: 7,032
作者:尹绍锴^*,闫理坦：东华大学理学院，上海
关键词:贝叶斯推断；分位数回归；有序数据；非对称拉普拉斯分布；吉布斯抽样；Bayesian Inference；Quantile Regression；Ordinal Data；Asymmetric Laplace Distribution；Gibbs Sampling

文章引用：尹绍锴, 闫理坦. 有序数据的贝叶斯分位数回归[J]. 统计学与应用, 2017, 6(5): 565-575. https://doi.org/10.12677/SA.2017.65064

1. 引言

1978年，Koenker and Bassett [1] 首次提出了用分位数回归方法来描述因变量的条件分位数与自变量之间的关系。分位数回归的提出引起了广泛的关注。分位数回归被广泛应用于农业、基因芯片技术、生存研究、经济学、医疗卫生、环境科学等领域。

在分位数回归估计中，算法的实现起着至关重要的作用。分位数回归模型估计方法中较经典的方法有Yu和Moyeed在2001年提出的单纯形法和内点法。其他的估计方法还有参数化方法和非参数方法。

有序数据在对事物分类的同时给出了各类的顺序，其数据仍表现为“类别”，但各类之间是有序的，可以比较优劣但其数值的大小差值却无意义。

在大多数研究领域中都会用到有序数据。目前已有大量的文献对有序数据进行了研究，但是用贝叶斯方法对有序数据进行分析目前还鲜少有人提及。

2. 贝叶斯分位数回归

对于一个随机变量x，如果它服从参数为 $\left(\mu ,\sigma ,p\right)$的非对称拉普拉斯分布(Asymmetric Laplace distribution, ALD)，那么其密度函数可以写为

$f\left(x;\mu ,\sigma ,p\right)=\frac{p\left(1-p\right)}{p}\exp \left\{-\frac{x-\mu }{\sigma }\left(p-1_{\left\{x\le \mu \right\}}\right)\right\},x\in R$(1)

其中 $$为位置参数， $\sigma >0$为尺度参数， $0<p<1$为偏度参数。

引理 如果随机变量 $\xi \sim ALD\left(\mu ,\sigma ,p\right)$则其方差 $\text{Var}\left(\xi \right)=\Psi \left(\sigma ,p\right)\ge 8{\sigma }^2$。

证明 由于

$E\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }xf\left(x;\mu ,\sigma ,p\right)\text{d}x}=\mu +\sigma \frac{1-2p}{p\left(1-p\right)}$

且

$\Psi \left(\sigma ,p\right)={\displaystyle {\int }_{-\infty }^{\infty }{\left(x-E\left(\xi \right)\right)}^{2}f\left(x;\mu ,\sigma ,p\right)\text{d}x}={\sigma }^2\frac{1-2p+2p^2}{p^2{\left(1-p\right)}^2}.$

根据函数

$g\left(p\right)=\frac{1-2p+2p^2}{p^2{\left(1-p\right)}^2}=\frac{1}{p^2}+\frac{1}{{\left(1-p\right)}^2}$

其中 $p\in \left(0,1\right)$，经计算可得

$\underset{p\in \left[0,1\right]}{\min }g\left(p\right)=g\left(\frac{1}{2}\right)=8$证毕。

如果随机变量 $\xi \sim ALD\left(\mu ,\sigma ,p\right)$，则其分位数函数(分布函数的逆函数)为

$F^{-1}\left(x;\mu ,\sigma ,p\right)=\{\begin{array}{ll}\mu +\frac{\sigma }{1-p}\log \frac{x}{p},\hfill & \text{if}0\le x\le p\hfill \\ \mu -\frac{\sigma }{p}\log \frac{1-x}{1-p},\hfill & \text{if}p<x\le 1\hfill \end{array}$(2)

从上式中可以的到一个很重要的性质：

${F^{-1}\left(x;\mu ,\sigma ,p\right)|}_{x=p}=\mu$(3)

即服从非对称拉普拉斯分布的随机变量在概率 $p$处的分位数等于位置参数 $\mu$。

设 $$服从如下的线性回归模型，

$r_t=Q_p\left(r_t|x_t,\beta \right)+u_t,t=1,2,\cdots ,n$(4)

式中 $Q_p\left(r_t|x_t,\beta \right)$是在观测到 $x_t$的条件下， $r_t$在 $p$概率水平下的分位数， $\beta$是参数向量，误差项 $\left\{u_t,t=1,2,\cdots ,n\right\}$相互独立且服从非对称拉普拉斯分布，即 $u_t\sim ALD\left(0,\sigma ,p\right)$。依据非对称拉普拉斯分布的

线性变换性质，式(4)可以等价表示为

$r_t~ALD\left(Q_p\left(r_t|x_t,\beta \right),\sigma ,p\right)$

故 $r_t$的密度函数为：

$f\left(r_t;Q_p\left(r_t|x_t,\beta \right),\sigma ,p\right)=\frac{p\left(1-p\right)}{p}\exp \left\{-\frac{Q_p\left(r_t|x_t,\beta \right)}{\sigma }\left[p-I\left(r_t\le Q_p\left(r_t|x_t,\beta \right)\right)\right]\right\}$(5)

$t=1,2,\cdots ,n$。样本的似然函数可以表示为

$L\left(r|\beta ,\sigma \right)=\frac{p^n{\left(1-p\right)}^n}{{\sigma }^n}\exp \left[-\frac{1}{\sigma }{\displaystyle {\sum }_{t=1}^n{\rho }_p\left(r_t-Q_p\left(r_t|x_t,\beta \right)\right)}\right]$(6)

其中 ${\rho }_p$为损失函数

${\rho }_p\left(u\right)=u\left(p-I\left(u<0\right)\right)>0$.

设参数 $\beta$和 $\sigma$相互独立，其先验密度分别为 $f\left(\beta \right)$和 $f\left(\sigma \right)$。根据贝叶斯定理，参数 $\beta ,\sigma$的联合后验密度为

$\pi \left(\beta ,\sigma |r\right)\propto L\left(r|\beta ,\sigma \right)f\left(\beta \right)f\left(\sigma \right)$(7)

在已知的文献中，如Yu等 [2] ，在估计过程中将尺度参数设定为常数1。但在实际应用中，将尺度参数设为常数缺有欠妥当。若尺度参数定为常数1，则随机变量的方差将不小8。在实际的数据研究中这个约束条件的存在并不合理。事实上，尺度参数的存在使得服从非对称拉普拉斯分布的随机变量的方差能取任意正值。尺度参数影响了参数估计的质量，在实际的应用当中，它应该被当为待估参数去处理，而不是主观地被设为某个常数。

对于分位数回归中的参数，它们不存在标准的共轭先验分布。这导致在贝叶斯推理中获取参数的后验分布解析表达式有一定的困难。若选择的参数的先验分布能确保参数的后验分布为真实分布，我们就可以利用MCMC 模拟得到参数的估计。

在Yu的研究中 [2] ，尺度参数为常数时，参数 $\beta$的先验分布是非真实均匀分布，得到的联合后验分布为真实分布。考虑到尺度参数的参数化，在无信息的前提下，参数的先验分布存在很多种可能，而被选择的先验分布可能会是非真实的分布。在这种情况下，参数 $\sigma$与 $\beta$满足什么条件会使得后验联合分布为真实分布。这对贝叶斯推理是很重要的一环。而下面的定理就可以保证后验分布是真实分布。

定理 对于任意 $\sigma >0$，参数的似然函数由式(6)给出， $\beta$和 $\sigma$的先验密度分别为 $f\left(\beta \right)$和 $f\left(\sigma \right)$。若 $\beta$服从不真实均匀分布， $f\left(\sigma \right)$满足

${\displaystyle {\int }_0^{\infty }f\left(\sigma \right)\frac{1}{{\sigma }^n}\text{d}\sigma }<\infty$

则参数的联合后验分布 $\pi \left(\beta ,\sigma |r\right)$是一个真实分布，即满足条件

$0<{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_0^{\infty }\pi \left(\beta ,\sigma |r\right)\text{d}\beta \text{d}\sigma }}<\infty$

证明 根据式(6)与(7)有

$\begin{array}{c}\pi \left(\beta ,\sigma |r\right)\propto L\left(r|\beta ,\sigma \right)f\left(\beta \right)f\left(\sigma \right)\\ =f\left(\beta \right)f\left(\sigma \right)\frac{p^n{\left(1-p\right)}^n}{{\sigma }^n}\left[-\frac{1}{\sigma }\underset{t=1}{\overset{n}{{\displaystyle \sum }}}{\rho }_p\left(r_t-Q_p\left(r_t|x_t,\beta \right)\right)\right]\end{array}$

可以推知

$0<\exp \left[-\frac{1}{\sigma }\underset{t=1}{\overset{n}{{\displaystyle \sum }}}{\rho }_p\left(r_t-Q_p\left(r_t|x_t,\beta \right)\right)\right]<1,0<{\displaystyle {\int }_0^{\infty }f\left(\sigma \right)\frac{1}{{\sigma }^n}\text{d}\sigma }<\infty ,$

于是

$0<L\left(r|\beta ,\sigma \right)f\left(\beta \right)f\left(\sigma \right)<f\left(\beta \right)f\left(\sigma \right)\frac{p^n{\left(1-p\right)}^n}{{\sigma }^n}$

非真实先验分布是在贝叶斯推断中在无信息先验情况下常采用的一种先验分布形式，参数 $\beta$服从非真实先验分布，则其概率密度的积分是无穷大的 [3] ，即

${\displaystyle {\int }_{-\infty }^{\infty }f\left(\beta \right)\text{d}\beta }=\infty$

可得

$\begin{array}{c}0={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_0^{\infty }0\text{d}\beta \text{d}\sigma }}\\ <{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_0^{\infty }L\left(r|\beta ,\sigma \right)f\left(\beta \right)f\left(\sigma \right)\text{d}\beta \text{d}\sigma }}\\ =\underset{i=1}{\overset{n}{{\displaystyle \prod }}}\underset{k=1}{\overset{K}{{\displaystyle \prod }}}{\left[F_{AL}\left(\frac{{\gamma }_{p,j}-{x^{\prime }}_i{\beta }_p}{\sigma }\right)-F_{AL}\left(\frac{{\gamma }_{p,j-1}-{x^{\prime }}_i{\beta }_p}{\sigma }\right)\right]}^{I\left(y_i=k\right)}\\ =p^n{\left(1-p\right)}^n{\displaystyle {\int }_{-\infty }^{\infty }f\left(\beta \right)\text{d}\beta }{\displaystyle {\int }_0^{\infty }f\left(\sigma \right)\frac{1}{{\sigma }^n}\text{d}\sigma }=\infty \end{array}$

因此 $0<{\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{\infty }\pi \left(\beta ,\sigma |r\right)\text{d}\beta \text{d}\sigma }<\infty$证毕。

3. 有序数据模型

有序数据的分位数回归模型能够用连续的潜变量 $t_i$来表示

$t_i={x^{\prime }}_i{\beta }_p+{ϵ}_i,i=1,\cdots ,n$(8)

其中， $x_i$是 $j$维的协向量， ${\beta }_p$是 $j$维的未知参数向量，误差项 ${ϵ}_i\sim ALD\left(0,1,p\right)$。 $n$表示观测值的个数。潜变量与我们观测到的响应变量 $y_i$存在一定的关系，且响应变量有 $K$个指标。我们用割点向量 ${\gamma }_p$来反应 $t_i$和 $y_i$的关系：

${\gamma }_{p,k-1}<t_i<{\gamma }_{p,k}\Rightarrow y_i=k,i=1,\cdots ,n;k=1,\cdots ,K$

其中 ${\gamma }_{p,0}=-\infty$， ${\gamma }_{p,K}=\infty$。一般而言 ${\gamma }_{p,1}$常被设定为0 [4] 。对于给定的数据 $y={\left(y_1,\cdots ,y_n\right)}^{\prime }$，模型的似然函数能够用包含有未知参数 $\left({\beta }_p,{\gamma }_p,\sigma \right)$的函数表示

$\begin{array}{c}f\left({\beta }_p,{\gamma }_p,\sigma ;y\right)=\underset{i=1}{\overset{n}{{\displaystyle \prod }}}\underset{k=1}{\overset{K}{{\displaystyle \prod }}}P{\left(y_i=k|{\beta }_p,{\gamma }_p,\sigma \right)}^{I\left(y_i=k\right)}\\ =\underset{i=1}{\overset{n}{{\displaystyle \prod }}}\underset{k=1}{\overset{K}{{\displaystyle \prod }}}{\left[F_{AL}\left(\frac{{\gamma }_{p,j}-{x^{\prime }}_i{\beta }_p}{\sigma }\right)-F_{AL}\left(\frac{{\gamma }_{p,j-1}-{x^{\prime }}_i{\beta }_p}{\sigma }\right)\right]}^{I\left(y_i=k\right)}\end{array}$

其中 $F_{AL}(\cdot )$表示非对称拉普拉斯分布的分布函数。 $I\left(y_i=k\right)$为示性函数。

为了在分位数回归中得到吉布斯算法，关于非对称拉普拉斯分布，Kozumi和Kobayashi在2011年 [5] 给出了一个基于指数分布和正态分布的混合表达式。若 ${ϵ}_i\sim ALD\left(0,1,p\right)$，则有

${ϵ}_i=\theta {\omega }_i+\tau \sqrt{{\omega }_i}{u}_i,i=1,\cdots ,n$(9)

其中 ${\omega }_i$和 $u_i$相互独立， $u_i\sim N\left(0,1\right)$， ${\omega }_i\sim E\left(1\right)$， $E$表示指数分布。式(9)中的常数项 $\left(\theta ,\tau \right)$为

$\theta =\frac{1-2p}{p\left(1-p\right)},\tau =\sqrt{\frac{2}{p\left(1-p\right)}}.$

非对称拉普拉斯分布的混合表达式使得在研究过程中，可以直接利用正态分布的性质。而这将对有序数据分位数回归模型的分析起到重大的作用。

贝叶斯方法来估计分位数回归模型运用到了潜变量表达式 (8)与非对称拉普拉斯分布的混合表达式(9)， $p$分位回归模型能够表示为

$t_i={x^{\prime }}_i{\beta }_p+\theta {\omega }_i+\tau \sqrt{{\omega }_i}{u}_i,i=1,\cdots ,n$(10)

上式中可以看出潜变量的条件分布服从正态分布，即 $t_i|{\beta }_p,{\omega }_i\sim N\left({x^{\prime }}_i{\beta }_p+\theta {\omega }_i,{\tau }^2{\omega }_i\right)$。这使得我们在估计

过程当中，能够利用正态分布的性质来进行研究。

考虑结果拥有三个指标的情况，回归模型能够用下式表示

$\{\begin{array}{l}{t}_i={x^{\prime }}_i{\beta }_p+{\sigma }_p{ϵ}_i={x^{\prime }}_i{\beta }_p+\theta {\omega }_i+\tau \sqrt{{\omega }_i}{u}_i,i=1,\cdots ,n\hfill \\ {\gamma }_{k-1}<t_i\le {\gamma }_k{y}_i=k,i=1,\cdots ,n\hfill \end{array}$(11)

其中 ${\sigma }_p$是 $p$分位水平下的尺度参数，割点向量 $\left({\gamma }_0,{\gamma }_1,{\gamma }_2,{\gamma }_3\right)$满足条件 ${\gamma }_0=-\infty ,{\gamma }_3=\infty$，且 ${\gamma }_1,{\gamma }_2$为某固定常值。

对于回归模型(11)，潜变量 $t_i$的条件期望包含了尺度参数 ${\sigma }_p$，这导致模型不能直接利用吉布斯抽样方法进行估计，但对模型进行简单的变形之后就能消除尺度参数的影响从而能利用吉布斯抽样方法来进行估计。模型(11) 能表示为

$t_i={x^{\prime }}_i{\beta }_p+\theta {\nu }_i+\tau \sqrt{{\sigma }_p{\nu }_i}{u}_i$(12)

其中 ${\nu }_i={\sigma }_p{\omega }_i$，因此潜变量的条件分布服从正态分布，即 $t_i|{\beta }_p,{\sigma }_p,{\nu }_i\sim N\left({x^{\prime }}_i{\beta }_p+\theta {\nu }_i,{\tau }^2{\sigma }_p{\nu }_i\right)$。接下来，

就是要设定参数的先验分布来获得参数的后验分布。设定参数的先验分布为：

$\{\begin{array}{l}{\beta }_p\sim N\left({\beta }_{p0},B_{p0}\right)\\ {\sigma }_p\sim IG\left(n_0/2,d_0/2\right)\\ {\nu }_i\sim E(\; \sigma \; p\; )\end{array}$

其中 $IG$和 $E$分别表示逆伽玛分布与指数分布。通过贝叶斯定理，参数 $\left(t,{\beta }_p,\nu ,{\sigma }_p\right)$的联合后验密度为：

$\begin{array}{c}\pi \left(t,{\beta }_p,\nu ,{\sigma }_p|y\right)\propto f\left(y|t,{\beta }_p,\nu ,{\sigma }_p\right)\pi \left(t|{\beta }_p,\nu ,{\sigma }_p\right)\pi \left(\nu |{\sigma }_p\right)\pi \left({\beta }_p\right)\pi \left({\sigma }_p\right)\\ \propto \left\{\underset{i=1}{\overset{n}{{\displaystyle \prod }}}f\left(y_i|t_i,{\sigma }_p\right)\right\}\pi \left(t|{\beta }_p,\nu ,{\sigma }_p\right)\pi \left(\nu |{\sigma }_p\right)\pi \left({\beta }_p\right)\pi \left({\sigma }_p\right).\end{array}$

当切割点与潜变量已知的情况下，观测值 $y_i$不依赖于 $\left({\beta }_p,\nu \right)$，根据式(12)，潜变量 $t_i$的条件密度 $\pi \left(t|{\beta }_p,\nu ,{\sigma }_p\right)$满足 $\pi \left(t|{\beta }_p,\nu ,{\sigma }_p\right)={\displaystyle {\prod }_{i=1}^{n}N\left({x^{\prime }}_i{\beta }_p+\theta {\nu }_i,{\tau }^2{\sigma }_p{\nu }_i\right)}$。待估参数的后验密度为

$\begin{array}{c}\pi \left(t,{\beta }_p,\nu ,{\sigma }_p|y\right)\propto \left\{\underset{i=1}{\overset{n}{{\displaystyle \prod }}}\underset{k=1}{\overset{3}{{\displaystyle \prod }}}I\left({\gamma }_{k-1}<t_i\le {\gamma }_k\right)N\left(t_i|{x^{\prime }}_i{\beta }_p+\theta {\nu }_i,{\tau }^2{\sigma }_p{\nu }_i\right)E\left({\nu }_i|{\sigma }_p\right)\right\}\\ \times N\left({\beta }_{p0},B_{p0}\right)IG\left(n_0/2,d_0/2\right),\end{array}$

据此，我们就能根据后验密度来导出我们感兴趣的参数的条件后验密度。

${\beta }_p$的条件后验密度 $\pi \left({\beta }_p|t,{\sigma }_p,\nu \right)$正比于 $\pi \left({\beta }_p\right)\times f\left(t|{\beta }_p,{\sigma }_p,\nu \right)$

$\begin{array}{c}\pi \left({\beta }_p|t,{\sigma }_p,\nu \right)\propto \exp \left[-\frac{1}{2}\left\{\underset{i=1}{\overset{n}{{\displaystyle \sum }}}{\left(\frac{t_i-{x^{\prime }}_i{\beta }_p-\theta {\nu }_i}{\tau \sqrt{{\sigma }_p{\nu }_i}}\right)}^2+{\left({\beta }_p-{\beta }_{p0}\right)}^{\prime }{B}_{p0}^{-1}\left({\beta }_p-{\beta }_{p0}\right)\right\}\right]\\ =\exp [-\frac{1}{2}\{{{\beta }^{\prime }}_p\left(\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\frac{x_i{x^{\prime }}_i}{{\tau }^2{\sigma }_p{\nu }_i}\right){\beta }_p-{{\beta }^{\prime }}_p\left(\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\frac{x_i\left(t_i-\theta {\nu }_i\right)}{{\tau }^2{\sigma }_p{\nu }_i}\right)-\left(\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\frac{{x^{\prime }}_i\left(t_i-\theta {\nu }_i\right)}{{\tau }^2{\sigma }_p{\nu }_i}\right)\\ +\left(\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\frac{{\left(t_i-\theta {\nu }_i\right)}^2}{{\tau }^2{\sigma }_p{\nu }_i}\right)+{\left({\beta }_p-{\beta }_{p0}\right)}^{\prime }{B}_{p0}^{-1}\left({\beta }_p-{\beta }_{p0}\right)\}]\\ \propto \exp [-\frac{1}{2}\{{{\beta }^{\prime }}_p\left(\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\frac{x_i{x^{\prime }}_i}{{\tau }^2{\sigma }_p{\nu }_i}+B_{p0}^{-1}\right){\beta }_p-{{\beta }^{\prime }}_p\left(\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\frac{x_i\left(t_i-\theta {\nu }_i\right)}{{\tau }^2{\sigma }_p{\nu }_i}+B_{p0}^{-1}\right)\\ -\left(\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\frac{{x^{\prime }}_i\left(t_i-\theta {\nu }_i\right)}{{\tau }^2{\sigma }_p{\nu }_i}+{{\beta }^{\prime }}_p{B}_{p0}^{-1}\right){\beta }_p\}]\end{array}$

设 ${\tilde{B}}_p^{-1}=\left({{\displaystyle \sum }}_{i=1}^n\frac{x_i{x^{\prime }}_i}{{\tau }^2{\sigma }_p{\nu }_i}+B_{p0}^{-1}\right)$， ${\tilde{\beta }}_p={\tilde{B}}_p\left({{\displaystyle \sum }}_{i=1}^n\frac{x_i\left(t_i-\theta {\nu }_i\right)}{{\tau }^2{\sigma }_p{\nu }_i}+B_{p0}^{-1}{\beta }_{p0}\right)$则有

$\pi \left({\beta }_p|t,{\sigma }_p,\nu \right)\propto \exp \left[-\frac{1}{2}\left\{{{\beta }^{\prime }}_p{\tilde{B}}_p^{-1}{\beta }_p-{{\beta }^{\prime }}_p{\tilde{B}}_p^{-1}{\tilde{\beta }}_p-{\widetilde{{\beta }^{\prime }}}_p{\tilde{B}}_p^{-1}{\beta }_p\right\}\right]$

经过简单的变形有

$\pi \left({\beta }_p|t,{\sigma }_p,\nu \right)\propto \exp \left[-\frac{1}{2}\left\{{{\beta }^{\prime }}_p{\tilde{B}}_p^{-1}{\beta }_p-{{\beta }^{\prime }}_p{\tilde{B}}_p^{-1}{\tilde{\beta }}_p-{{\tilde{\beta }}^{\prime }}_p{\tilde{B}}_p^{-1}{\beta }_p+{{\tilde{\beta }}^{\prime }}_p{\tilde{B}}_p^{-1}{\tilde{\beta }}_p-{{\tilde{\beta }}^{\prime }}_p{\tilde{B}}_p^{-1}{\tilde{\beta }}_p\right\}\right]$

因此

$\pi \left({\beta }_p|t,{\sigma }_p,\nu \right)\propto \exp \left[-\frac{1}{2}\left\{{\left({\beta }_p-{\tilde{\beta }}_p\right)}^{\prime }{\tilde{B}}_p^{-1}\left({\beta }_p-{\tilde{\beta }}_p\right)\right\}\right]$

结果表明， ${\beta }_p$的条件后验分布为正态分布，即 ${\beta }_p|t,{\sigma }_p,\nu \sim N\left({\tilde{\beta }}_p,{\tilde{B}}_p\right)$。

尺度参数 ${\sigma }_p$的条件后验密度 $\pi \left({\sigma }_p|t,{\beta }_p,\nu \right)$正比于 $f\left(t|{\beta }_p,\nu ,{\sigma }_p\right)\pi \left(\nu |{\sigma }_p\right)\pi (\; \sigma \; p\; )$

$\begin{array}{c}\pi \left({\sigma }_p|t,{\beta }_p,\nu \right)\propto \underset{i=1}{\overset{n}{{\displaystyle \prod }}}\left\{{\sigma }_p^{-1/2}\exp \left[-\frac{1}{2}{\left(\frac{t_i-{x^{\prime }}_i{\beta }_p-\theta {\nu }_i}{\tau \sqrt{{\sigma }_p{\nu }_i}}\right)}^2\right]\times {\sigma }_p^{-1}\exp \left(-\frac{{\nu }_i}{{\sigma }_p}\right)\right\}\exp \left[-\frac{d_0}{2{\sigma }_p}\right]{\sigma }_p^{-\left(\frac{n_0}{2}+1\right)}\\ \propto {\sigma }_p^{-\left(\frac{n_0}{2}++\frac{3n}{2}+1\right)}\exp \left[-\frac{1}{{\sigma }_p}\left\{\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\frac{{\left(t_i-{x^{\prime }}_i{\beta }_p-\theta {\nu }_i\right)}^2}{2{\tau }^2{\nu }_i}+\frac{d_0}{2}+\underset{i=1}{\overset{n}{{\displaystyle \sum }}}{\nu }_i\right\}\right]\end{array}$

其中 $\tilde{n}=n_0+3n$， $\tilde{d}={{\displaystyle \sum }}_{i=1}^n{\left(t_i-{x^{\prime }}_i\beta -\theta {\nu }_i\right)}^2/{\tau }^2{\nu }_i+d_0+2{{\displaystyle \sum }}_{i=1}^n{\nu }_i$。 ${\sigma }_p$的条件后验分布为逆伽玛分布即 ${\sigma }_p|t,{\beta }_p,\nu \sim IG\left(\tilde{n}/2,\tilde{d}/2\right)$。

$\nu$的条件后验密度 $\pi \left(\nu |t,{\beta }_p,{\sigma }_p\right)$正比于 $f\left(t|{\beta }_p,\nu ,{\sigma }_p\right)\pi \left(\nu \right)$且 ${\nu }_i$有如下形式

$\begin{array}{c}\pi \left({\nu }_i|w,{\beta }_p,{\sigma }_p\right)\propto {\nu }_i^{-1/2}\exp \left[-\frac{1}{2}{\left(\frac{w_i-{x^{\prime }}_i{\beta }_p-\theta {\nu }_i}{\tau \sqrt{{\sigma }_p{\nu }_i}}\right)}^2-\frac{{\nu }_i}{{\sigma }_p}\right]\\ \propto {\nu }_i^{-1/2}\exp \left[-\frac{1}{2{\sigma }_p}\left(\frac{{\left(w_i-{x^{\prime }}_i{\beta }_p\right)}^2+{\theta }^2{\nu }_i^2-2\theta {\nu }_i\left(t_i-{x^{\prime }}_i{\beta }_p\right)}{{\tau }^2{\nu }_i}+2{\nu }_i\right)\right]\\ ={\nu }_i^{-1/2}\exp \left[-\frac{1}{2}\left\{\frac{{\left(t_i-{x^{\prime }}_i{\beta }_p\right)}^2}{{\tau }^2{\sigma }_p}{\nu }_i^{-1}+\left(\frac{{\theta }^2}{{\tau }^2{\sigma }_p}+\frac{2}{{\sigma }_p}\right){\nu }_i-\frac{2\theta \left(t_i-{x^{\prime }}_i{\beta }_p\right)}{{\tau }^2{\sigma }_p}\right\}\right]\\ \propto {\nu }_i^{-1/2}\exp \left[-\frac{1}{2}\left\{\frac{{\left(t_i-{x^{\prime }}_i{\beta }_p\right)}^2}{{\tau }^2{\sigma }_p}{\nu }_i^{-1}+\left(\frac{{\theta }^2}{{\tau }^2{\sigma }_p}+\frac{2}{{\sigma }_p}\right){\nu }_i\right\}\right]\end{array}$

设 ${\tilde{\lambda }}_i=\frac{{\left(t_i-{x^{\prime }}_i{\beta }_p\right)}^2}{{\tau }^2{\sigma }_p}$， $\tilde{\eta }=\left(\frac{{\theta }^2}{{\tau }^2{\sigma }_p}+\frac{2}{{\sigma }_p}\right)$则有

$\pi \left({\nu }_i|t,{\beta }_p,{\sigma }_p\right)\propto {\nu }_i^{-1/2}\exp \left[-\frac{1}{2}\left\{{\tilde{\lambda }}_i{\nu }_i^{-1}+\tilde{\eta }{\nu }_i\right\}\right]$

${\nu }_i$的条件后验分布为广义逆伽玛分布，即 ${\nu }_i|t,{\beta }_p,{\sigma }_p\sim GIG\left(0.5,{\tilde{\lambda }}_i,\tilde{\eta }\right)$。

4. 数据模拟分析

我们从如下模型获得观测值

$\{\begin{array}{l}{Y}_i={\beta }_0+{\beta }_{1}x+\nu ,i=1,2,\cdots ,n,\hfill \\ x\sim unif\left(0,10\right),\hfill \\ {\beta }_0=1,{\beta }_1=2.\hfill \end{array}$

${\beta }_0$与 ${\beta }_1$的先验分布为 $N\left(0,10\right)$，尺度参数 $\sigma$的先验分布设定自由度为3的卡方分布。采用Gibbs抽样法共模拟5000次，为消除初值对抽样分布的影响，去掉前2000个抽样值。在应用Gibbs抽样算法对参数进行抽样时，特定参数的Gibbs模拟值是否是真实后验分布的合理近似，这对参数估计值的正确性意义重大，因此在统计推断之前需要对抽样分布进行检验。

图1给出了样本容量为100时， ${\beta }_0$， ${\beta }_1$和 $\sigma$在0.5分位点下的MCMC抽样值轨迹和核密度图。图中显示三个参数的3000个抽样值序列在某一个值的上下波动，因此抽样构成的马尔科夫链收敛

图2给出了样本容量为100时， ${\beta }_0$， ${\beta }_1$和 $\sigma$在0.5分位点下抽样值的自相关图。从图中可以看出随着滞后期的增加，自相关系数趋向于零。

综上，马尔科夫链收敛且自相关系数趋于零，所以在样本容量为100，各参数的贝叶斯中位数回归估计结果是可靠的。其他样本容量和分位点下的马尔科夫链均收敛，自相关系数也趋于零。表1和表2给出了不同样本容量下得到的MCMC参数估计结果。

在表1和表2中。比较 ${\beta }_0$， ${\beta }_1$和 $\sigma$参数化与未参数化下的估计结果，可以发现在一定的样本容量，特定的分位点下，当 $\sigma$被参数化时， ${\beta }_0$和 ${\beta }_1$的抽样标准差均小于 $\sigma$未参数化是各自的标准差。这表明尺度参数的参数化可以提高被估系数的精确度，从而提高待估系数假设检验的准确性。因此在用Gibbs抽样方法进行贝叶斯分位数回归估计是应该将尺度参数参数化，而不是设定为常数。此外，还可以发现，

表1. ${\beta }_0$在不同分位数水平下的估计量

表2. ${\beta }_1$在不同分位数水平下的估计量

注：给定样本容量，未参数化和参数化情形下的数字特征中，第一行表示参数的真值，第二行和第三行数值分别表示参数估计量分布的均值和标准差。

在给定样本容量的情况下，无论 $\sigma$是否参数化， ${\beta }_1$的标准差总是小于同一分位数水平下 ${\beta }_0$的标准差。对于 ${\beta }_1$，无论尺度参数是否被参数化，它的精度都随着样本量的增加而提高。

参考文献

[1]	Koenker, R. and Bassett, G. (1978) Regression Quantiles. Econometrica, 46, 33-50. https://doi.org/10.2307/1913643
[2]	Yu, K. and Moyeed, R.A. (2001) Bayesian Quantile Regression. Statistics and Probability Letters, 54, 437-447. https://doi.org/10.1016/S0167-7152(01)00124-9
[3]	Zellner, A. (1971) An Introduction to Bayesian Inference in Econometrics. John Wiley Press, New York.
[4]	Jeliazkov, I., Graves, J. and Kutzbach, M. (2008). Fitting and Comparison of Models for Multivariate Ordinal Out- comes. Advances in Econometrics: Bayesian Econometrics, 23, 115-156. https://doi.org/10.1016/S0731-9053(08)23004-5
[5]	Kozumi, H. and Kobayashi, G. (2011) Gibbs Sampling Methods for Bayesian Quantile Regression. Journal of Statistical Computation and Simulation, 81, 1565-1578. https://doi.org/10.1080/00949655.2010.496117