贝叶斯线性模型
考虑多元正态线性回归模型
$y_i=\beta_0+\beta_1x_{i1}+\cdots+\beta_{p-1}x_{i,p-1}+\varepsilon_i,\quad\varepsilon_i\sim N(0,\sigma^2).$
记 $\pmb y=(y_1,\dots,y_n)’$，$\pmb X$ 为 $n\times p$ 设计矩阵，$\pmb\beta=(\beta_0,\dots,\beta_{p-1})’$ 为感兴趣的参数，则参数的最小二乘估计（也是极大似然估计）
$\hat{\pmb\beta}=(\pmb X'\pmb X)^{-1}\pmb X'\pmb y\sim N(\pmb\beta,\sigma^2(\pmb X'\pmb X)^{-1}).$
对方差 $\sigma^2$ 的估计为
$\hat{\sigma}^2=\dfrac{1}{n-p}(\pmb y-\pmb X\hat{\pmb\beta})'(\pmb y-\pmb X\hat{\pmb\beta}).$
后验分布
在Lec2.3节中我们已经详细推导了无信息先验与共轭先验下的参数后验分布，这里我们对不正常先验和NIG先验作简要的回顾（另外，Lec8.2节已经提供了Zellner’s g先验的结果）。
由模型可知 $(\pmb\beta,\sigma^2)$ 的似然函数为
$\begin{array}{rl}
L(\pmb\beta,\sigma^2\mid\pmb y,\pmb X)\propto&\!\!\!\!(\sigma^2)^{-\frac{n}{2}}\exp\left\{-\dfrac{(\pmb y-\pmb X\pmb\beta)'(\pmb y-\pmb X\pmb\beta)}{2\sigma^2}\right\}\\
\propto&\!\!\!\!(\sigma^2)^{-\frac{n}{2}}\exp\left\{-\dfrac{(\pmb\beta-\hat{\pmb\beta})'\pmb X'\pmb X(\pmb\beta-\hat{\pmb\beta})+\nu s^2}{2\sigma^2}\right\},
\end{array}$
其中 $\nu=n-p$，$\nu s^2=(\pmb y-\pmb X\hat{\pmb\beta})’(\pmb y-\pmb X\hat{\pmb\beta})$。
对 $\pmb\beta$ 采用Jeffreys先验，即 $\pi(\pmb\beta)=1$，此时若 $\sigma^2$ 已知，则
$\pmb\beta\mid\pmb y,\pmb X\sim N(\hat{\pmb\beta},\sigma^2(\pmb X'\pmb X)^{-1}).$
大部分情况下 $\sigma^2$ 是未知的，此时若 $\pi(\sigma^2)=1/\sigma^2$，且 $\pmb\beta$ 与 $\sigma^2$ 独立，则
$\pmb\beta\mid\pmb y,\pmb X\sim\mathscr{T}_p\left(\nu,\hat{\pmb\beta},s^2(\pmb X'\pmb X)^{-1}\right),\quad\sigma^2\mid\pmb y,\pmb X\sim\Gamma^{-1}\left(\dfrac{\nu}{2},\dfrac{\nu s^2}{2}\right).$
若假设 $\sigma^2$ 是共轭先验，即 $\sigma^2\sim\Gamma^{-1}(a,b)$，且 $\pmb\beta$ 与 $\sigma^2$ 独立，则
$\pmb\beta\mid\pmb y,\pmb X\sim\mathscr{T}_p\left(\nu+2a,\hat{\pmb\beta},\dfrac{\nu s^2+2b}{\nu+2a}(\pmb X'\pmb X)^{-1}\right),\quad\sigma^2\mid\pmb y,\pmb X\sim\Gamma^{-1}\left(\dfrac{\nu}{2}+a,\dfrac{\nu s^2}{2}+b\right).$
实际中，$a$ 和 $b$ 都取较小值（如0.001），相当于前一种情形。
设 $\pmb\beta\mid\sigma^2\sim N_p(\pmb\mu,\sigma^2\pmb V)$，$\sigma^2\sim\Gamma^{-1}(a,b)$（比Lec2.3中给的 $\pmb\beta$ 分量独立的情况更加一般），则
$\begin{array}{rl}
\pi(\pmb\beta,\sigma^2\mid\pmb y,\pmb X)\propto&\!\!\!\!(\sigma^2)^{-\frac{n}{2}}\exp\left\{-\dfrac{(\pmb\beta-\hat{\pmb\beta})'\pmb X'\pmb X(\pmb\beta-\hat{\pmb\beta})+\nu s^2}{2\sigma^2}\right\}\\
&\quad\cdot(\sigma^2)^{-\frac{p}{2}}\exp\left\{-\dfrac{(\pmb\beta-\pmb\mu)'\pmb V^{-1}(\pmb\beta-\pmb\mu)}{2\sigma^2}\right\}\cdot(\sigma^2)^{-a-1}\exp\left\{-\dfrac{b}{\sigma^2}\right\}\\
\propto&\!\!\!\!(\sigma^2)^{-\frac{n+p}{2}-a-1}\exp\left\{-\dfrac{(\pmb\beta-\tilde{\pmb\mu})'\tilde{\pmb V}^{-1}(\pmb\beta-\tilde{\pmb\mu})+\tilde{b}}{2\sigma^2}\right\},
\end{array}$
其中 $\tilde{\pmb\mu}=(\pmb X’\pmb X+\pmb V^{-1})^{-1}(\pmb X’\pmb y+\pmb V^{-1}\pmb\mu)$，$\tilde{\pmb V}=(\pmb X’\pmb X+\pmb V^{-1})^{-1}$，
$\begin{array}{rl}
\tilde{b}=&\!\!\!\!2b+\pmb y'\pmb y+\pmb\mu'\pmb V^{-1}\pmb\mu-\tilde{\pmb\mu}'\tilde{\pmb V}^{-1}\tilde{\pmb\mu}\\
=&\!\!\!\!2b+\pmb y'\pmb y+\pmb\mu'\pmb V^{-1}\pmb\mu-(\pmb X'\pmb y+\pmb V^{-1}\pmb\mu)'\tilde{\pmb V}(\pmb X'\pmb y+\pmb V^{-1}\pmb\mu)\\
=&\!\!\!\!2b+\pmb y'(\pmb I_n-\pmb X\tilde{\pmb V}\pmb X')\pmb y-2\pmb y'\pmb X\tilde{\pmb V}\pmb V^{-1}\pmb\mu+\pmb\mu'(\pmb V^{-1}-\pmb V^{-1}\tilde{\pmb V}\pmb V^{-1})\pmb\mu
\end{array}$
注意到恒等式 $(\pmb A+\pmb B\pmb D\pmb C)^{-1}=\pmb A^{-1}-\pmb A^{-1}\pmb B(\pmb D^{-1}+\pmb C\pmb A^{-1}\pmb B)^{-1}\pmb C\pmb A^{-1}$，于是
$\pmb V^{-1}-\pmb V^{-1}\tilde{\pmb V}\pmb V^{-1}=[(\pmb X'\pmb X)^{-1}+\pmb V]^{-1}=\pmb X'\pmb X-\pmb X'\pmb X\tilde{\pmb V}\pmb X'\pmb X=\pmb X'(\pmb I_n-\pmb X\tilde{\pmb V}\pmb X')\pmb X.$
又 $\tilde{\pmb V}(\pmb X’\pmb X+\pmb V^{-1})=\pmb I_p$，所以 $\pmb X\tilde{\pmb V}\pmb V^{-1}=(\pmb I_n-\pmb X\tilde{\pmb V}\pmb X’)\pmb X$，由上述结果可知
$\begin{array}{rl}
\tilde{b}=&\!\!\!\!2b+(\pmb y-\pmb X\pmb\mu)'(\pmb I_n-\pmb X\tilde{\pmb V}\pmb X')(\pmb y-\pmb X\pmb\mu)\\
=&\!\!\!\!2b+(\pmb y-\pmb X\pmb\mu)'(\pmb I_n+\pmb X\pmb V\pmb X')^{-1}(\pmb y-\pmb X\pmb\mu).
\end{array}$
这里第二个等号又用到了前面的矩阵恒等式。所以
$\pmb\beta\mid\pmb y,\pmb X\sim\mathscr{T}_p\left(n+2a,\tilde{\pmb\mu},\dfrac{\tilde{b}}{n+2a}\tilde{\pmb V}\right),\quad\sigma^2\mid\pmb y,\pmb X\sim\Gamma^{-1}\left(\dfrac{n}{2}+a,\dfrac{\tilde{b}}{2}\right).$
边际分布
注意到模型可以写为
$\pmb y=\pmb X\pmb\beta+\pmb\varepsilon_1,\quad\pmb\beta=\pmb\mu+\pmb\varepsilon_2,$
其中 $\pmb\varepsilon_1\sim N(\pmb 0,\sigma^2\pmb I)$，$\pmb\varepsilon_2\sim N(\pmb 0,\sigma^2\pmb V)$，而 $\sigma^2\sim\Gamma^{-1}(a,b)$。所以
$\pmb y=\pmb X\pmb\mu+\pmb X\pmb\varepsilon_2+\pmb\varepsilon_1\sim N(\pmb X\pmb\mu,\sigma^2(\pmb I+\pmb X\pmb V\pmb X')).$
因此
$\begin{array}{rl}
p(\pmb y)=&\!\!\!\!\displaystyle\int_{\Sigma^2}p(\pmb y\mid\sigma^2)\pi(\sigma^2)\text{d}(\sigma^2)\\
\propto&\!\!\!\!\displaystyle\int_{\Sigma^2}(\sigma^2)^{-\frac{n}{2}-a-1}\exp\left\{-\frac{(\pmb y-\pmb X\pmb\mu)'(\pmb I+\pmb X\pmb V\pmb X')^{-1}(\pmb y-\pmb X\pmb\mu)+2b}{2\sigma^2}\right\}\text{d}(\sigma^2)\\
\propto&\!\!\!\!\left[1+\dfrac{1}{2a}(\pmb y-\pmb X\pmb\mu)'\left(\dfrac{b}{a}(\pmb I+\pmb X\pmb V\pmb X')\right)^{-1}(\pmb y-\pmb X\pmb\mu)\right]^{-\frac{n+2a}{2}},
\end{array}$
于是 $\pmb y$ 的边际分布为
$\pmb y\sim\mathscr{T}_n\left(2a,\pmb X\pmb\mu,\dfrac{b}{a}(\pmb I+\pmb X\pmb V\pmb X')\right).$
预测分布
当有了一组新观测 $\pmb X_{\text{new}}$ 后，我们希望预测其响应 $\pmb y_{\text{new}}$。由预测公式
$\begin{array}{rl}
p(\pmb y_{\text{new}}\mid\pmb y)=&\!\!\!\!\displaystyle\iint p(\pmb y_{\text{new}}\mid\pmb\beta,\sigma^2)\pi(\pmb\beta,\sigma^2\mid\pmb y,\pmb X)\text{d}\pmb\beta\text{d}(\sigma^2)\\
=&\!\!\!\!\displaystyle\iint N(\pmb X_{\text{new}}\pmb\beta,\sigma^2\pmb I_m)\cdot NIG\left(\tilde{\pmb\mu},\tilde{\pmb V},\dfrac{n+2a}{2},\dfrac{\tilde{b}}{2}\right)\text{d}\pmb\beta\text{d}(\sigma^2)\\
=&\!\!\!\!\mathscr{T}_m\left(n+2a,\pmb X_{\text{new}}\tilde{\pmb\mu},\dfrac{\tilde{b}}{n+2a}(\pmb I+\pmb X_{\text{new}}\tilde{\pmb V}\pmb X_{\text{new}}')\right),
\end{array}$
注意这里利用了正态与NIG的边际分布结论。
特别地，当取无信息先验 $\pi(\pmb\beta,\sigma^2)=1/\sigma^2$ 时，等价于在NIG系统中令 $\pmb V^{-1}\rightarrow\pmb 0$，$a\rightarrow-p/2$，$b\rightarrow0$，此时，预测分布为
$\pmb y_{\text{new}}\sim\mathscr{T}_m\left(n-p,\pmb X_{\text{new}}\hat{\pmb\beta},s^2(\pmb I_m+\pmb X_{\text{new}}(\pmb X'\pmb X)^{-1}\pmb X_{\text{new}})\right).$
回归诊断
与频率统计下类似的，我们考虑残差 $\pmb e=\pmb y-\pmb X\pmb\beta$ 为一个我们希望推断的未知量。如果 $\pmb e$ 已知，则可以用来检查模型假设、独立性等。
记 $\hat{\pmb e}$ 为 $\pmb e$ 的后验均值，即
$\hat{\pmb e}=\mathbb E(\pmb e\mid\pmb y,\pmb X)=\pmb y-\pmb X\cdot\mathbb E(\pmb\beta\mid\pmb y,\pmb X).$
在无信息先验下，$\mathbb E(\pmb\beta\mid\pmb y,\pmb X)=\hat{\pmb\beta}=(\pmb X’\pmb X)^{-1}\pmb X’\pmb y$。根据后验分布的结论可知，
$t=\dfrac{e_i-\hat{e}_i}{s_{e_i}}\sim t(n-p),$
其中 $s_{e_i}^2=\pmb x_i’(\pmb X’\pmb X)^{-1}\pmb x_is^2$，$e_i=y_i-\pmb x_i’\pmb\beta$ 为 $\pmb e$ 的第 $i$ 个分量。进而 $e_i$ 的 $1-\alpha$ 贝叶斯可信区间为
$\left(\hat{e}_i-t_{\alpha/2}(n-p)\cdot s_{e_i},\hat{e}_i+t_{\alpha/2}(n-p)\cdot s_{e_i}\right).$
通过观察各个点是否在此可信区间里来判断模型是否合适。
对于影响点（或异常点）分析，可以通过对每个异常点引入或剔除时模型系数估计的变化来判断。更一般地，可以将数据分成不同的子集，研究每个子集下参数估计的变化情况。
贝叶斯广义线性回归
广义线性模型（GLM）是线性模型的一种推广，它由三个部分组成：
下表给出了一些分布下的连接函数及其逆。
| 分布 |
连接函数 $\theta=g(\mu)$ |
逆连接函数 $\mu=g^{-1}(\theta)$ |
| 泊松 |
$\log(\mu)$ |
$\exp(\theta)$ |
| 二项 |
logit连接：$\log\dfrac{\mu}{1-\mu}$ |
$\dfrac{\exp(\theta)}{1+\exp(\theta)}$ |
| |
probit连接：$\Phi^{-1}(\mu)$ |
$\Phi(\theta)$ |
| |
cloglog连接：$\log(-\log(1-\mu))$ |
$1-\exp(-\exp(\theta))$ |
| 正态 |
$\mu$ |
$\theta$ |
| 伽马 |
$-\dfrac{1}{\mu}$ |
$-\dfrac{1}{\theta}$ |
| 负二项 |
$\log(1-\mu)$ |
$1-\exp(\theta)$ |
Deviance Function
尽管残差可以表示成类似线性模型中观测结果减去预测结果值，具体有
- 线性模型残差：$\pmb R_{\text{standard}}=\pmb y-\pmb X\pmb\beta$；
- 响应变量残差：$\pmb R_{\text{response}}=\pmb y-g^{-1}(\pmb X\pmb\beta)$；
- Pearson残差：$\pmb R_{\text{Pearson}}=\dfrac{\pmb y-\pmb\mu}{\sqrt{\text{Var}(\pmb\mu)}}$；
- Working残差：$\pmb R_{\text{working}}=(\pmb y-\pmb\mu)\dfrac{\partial\pmb\mu}{\partial\pmb\theta}=(\pmb y-\pmb\mu)/g’(\pmb\mu)$。
但是在GLM中，更有用的是使用Deviance residuals：
- 个体偏差函数：$R_{\text{Deviance}}=\text{sign}(y_i-\mu_i)\sqrt{\vert d(\pmb\theta,y_i)\vert}$；
- 总偏差：$\displaystyle D=\sum_{i=1}^n d(\pmb\theta,y_i)$；
其中 $d(\pmb\theta,y_i)=-2\left[l(\hat{\pmb\theta},\phi\mid y_i)-l(\tilde{\pmb\theta},\phi\mid y_i)\right]$，$\tilde{\pmb\theta}$ 为饱和模型参数估计。
下表给出了一些分布下的典则参数与偏差函数。
| 分布 |
典则参数 |
偏差函数 |
| 泊松 $\mathcal P(\mu)$ |
$\theta=\log(\mu)$ |
$\displaystyle 2\sum_{i=1}^n\left[y_i\log\dfrac{y_i}{\mu_i}-y_i+\mu_i\right]$ |
| 二项 $B(m,p)$ |
$\theta=\log\dfrac{\mu}{1-\mu}$ |
$\displaystyle 2\sum_{i=1}^n\left[y_i\log\dfrac{y_i}{\mu_i}+(m_i-y_i)\log\dfrac{m_i-y_i}{m_i-\mu_i}\right]$ |
| 正态 $N(\mu,\sigma^2)$ |
$\theta=\mu$ |
$\displaystyle\sum_{i=1}^n(y_i-\mu_i)^2$ |
| 伽马 $\Gamma(\mu,\delta)$ |
$\theta=-\dfrac{1}{\mu}$ |
$\displaystyle 2\sum_{i=1}^n\left[-\dfrac{y_i-\mu_i}{\mu_i}\log\dfrac{y_i}{\mu_i}\right]$ |
| 负二项 $NB(\mu,p)$ |
$\theta=\log(1-\mu)$ |
$\displaystyle 2\sum_{i=1}^n\left[y_i\log\dfrac{y_i}{\mu_i}+(1+y_i)\log\dfrac{1+\mu_i}{1+y_i}\right]$ |
下面我们进行GLM的贝叶斯推断。注意后验分布为
$\pi(\pmb\beta,\phi)\propto f(\pmb y;\pmb X,\pmb\beta,\phi)\pi(\pmb\beta,\phi)=\exp\left\{\dfrac{\pmb y'\pmb\theta-b(\pmb\theta)}{a(\phi)}+c(\pmb y,\phi)\right\}\pi(\pmb\beta,\phi).$
这里先验可以取无信息的，也可以是有信息的，不过一般来说，后验分布并没有简单的形式。
设 $Y\sim B(n,\pi)$，其中 $\pi$ 为成功概率，考虑 $\theta=\alpha+\pmb x’\pmb\beta=g(\pi)$，常用的连接函数 $g$ 有
$g(\pi)=\text{logit}(\pi)=\log\dfrac{\pi}{1-\pi},$ 或 $g(\pi)=\Phi^{-1}(\pi).$
当给定一组样本（样本量为 $N$）后，似然函数为
$\displaystyle\prod_{i=1}^N\binom{n}{y_i}\pi^{y_i}(1-\pi)^{n-y_i}=\prod_{i=1}^N\binom{n}{y_i}\left(\dfrac{\text{e}^\theta}{1+\text{e}^\theta}\right)^{y_i}\left(\dfrac{1}{1+\text{e}^\theta}\right)^{n-y_i},$
或
$\displaystyle\prod_{i=1}^N\binom{n}{y_i}\pi^{y_i}(1-\pi)^{n-y_i}=\prod_{i=1}^N\binom{n}{y_i}\left(\Phi(\theta)\right)^{y_i}\left(1-\Phi(\theta)\right)^{n-y_i}.$
根据问题背景，可以使用各种有信息或无信息先验分布。在logit连接函数下，研究表明一个厚尾分布比较合适，包括 $t$ 分布、柯西分布等。
设 $Y\sim\mathcal P(\lambda)$，考虑 $\log(\lambda)=\pmb x’\pmb\beta$，进而可以实现贝叶斯推断。
考虑飞机受损数据（Montgomery et al.，2006），数据中包含了越南战争中30次空袭任务中飞机受损的数据，共有4个变量：
- damage：飞机上损坏位置的个数；
- type：0-1变量表示飞机的类型（0表示A4；1表示A6）；
- bombload：载弹情况（吨）；
- airexp：飞行员飞行经验总月数。
使用Poisson对数回归模型
$\begin{array}{rl}
\text{damage}_i\sim&\!\!\!\!\mathcal P(\lambda_i),\quad i=1,2,\dots,30,\\
\log\lambda_i=&\!\!\!\!\beta_1+\beta_2\text{type}_i+\beta_3\text{bombload}_i+\beta_4\text{airexp}_i,
\end{array}$
下表给出了估计结果。
| node |
mean |
sd |
MC error |
2.5% |
median |
97.5% |
harmonic |
| $\beta_1$ |
-0.766 |
1.089 |
0.1762 |
-3.168 |
-0.835 |
1.619 |
|
| $\beta_2$ |
0.580 |
0.466 |
0.0513 |
-0.302 |
0.584 |
1.537 |
|
| $\beta_3$ |
0.177 |
0.068 |
0.0099 |
0.040 |
0.177 |
0.308 |
|
| $\beta_4$ |
-0.011 |
0.010 |
0.0015 |
-0.033 |
-0.010 |
0.007 |
|
| $B_1$ |
0.862 |
1.221 |
0.1829 |
0.042 |
0.434 |
5.050 |
0.465 |
| $B_2$ |
1.993 |
0.996 |
0.1050 |
0.739 |
1.793 |
4.652 |
1.786 |
| $B_3$ |
1.197 |
0.081 |
0.0118 |
1.041 |
1.193 |
1.360 |
1.194 |
| $B_4$ |
0.989 |
0.010 |
0.0015 |
0.968 |
0.990 |
1.007 |
0.989 |
其中 $B_j$ 的调和均值（最后一列）使用 $\beta_j$ 的后验均值额外计算。
从结果中可以看到
- 只有bombload的系数是远离0的，事实上根据DIC选择最优模型也是仅包含bombload的模型最优。
- A6的期望受损数是A4的两倍（当飞行员飞行经验时长和飞机载弹量相同），相同飞行员经验时长和载弹量情况下，A6比A4的受损程度从后验上预期高79%。
- 每吨载弹增加20%飞机的期望受损位置数。
- 飞行经验多一年减少1%的飞机期望受损位置数。
贝叶斯惩罚回归
当自变量个数比较多时，惩罚回归技术被广泛用来防止过度拟合。惩罚回归能够从大量自变量中选出能够预测因变量的自变量，因此常用于高维回归问题中。
惩罚回归的基本想法是对目标损失函数加上惩罚项，以使得比较小的系数可以压缩为0，而保留大的系数。例如
$\underset{\beta_0,\pmb\beta}{\text{minimize}}~\dfrac{1}{2n}\|\pmb y-\beta_0\pmb 1-\pmb X\pmb\beta\|_2^2+\lambda\|\pmb\beta\|_q,$
其中 $|\pmb\beta|_q$ 是 $\pmb\beta$ 的 $q$-范数，$\lambda$ 为惩罚参数，其值越大导致压缩到0的程度更大。当 $q=1$ 时即为著名的Lasso回归。
在贝叶斯回归 $y_i\mid x_i,\beta_0,\pmb\beta,\sigma^2\sim N(\beta_0+\pmb x_i’\pmb\beta,\sigma^2)$ 中，关于 $\pmb\beta$ 的负对数后验密度为
$\displaystyle -\log p(\pmb\beta\mid\pmb y,\pmb X)=\dfrac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\beta_0-\pmb x_i'\pmb\beta)^2-\log p(\pmb\beta\mid\lambda).$
参数的先验分布可以选为
$p(\beta_0,\pmb\beta,\sigma^2,\lambda)=p(\beta_0)p(\pmb\beta\mid\sigma^2,\lambda)p(\sigma^2)p(\lambda).$
注意这里 $\pmb\beta$ 的先验分布依赖于 $\sigma^2$ 和 $\lambda$，因为在某些情形下为了保证后验是单峰的，其需要依赖于 $\sigma^2$。
此外 $\lambda$ 是超参数，其作用类似于经典惩罚回归中的惩罚参数作用。对 $\lambda$ 的选取，有几种做法：
下表给出了贝叶斯惩罚回归常见压缩先验。
Ridge
从形式上可以看到，岭回归的惩罚项相应于先验为
$\begin{array}{rl}
\beta_j\mid\sigma^2,\lambda\sim&\!\!\!\!N\left(0,\dfrac{\sigma^2}{\lambda}\right),\quad j=1,\dots,p,\\
\lambda\sim&\!\!\!\!C_+(0,1),
\end{array}$
这里 $\lambda$ 的先验选择了半柯西先验 $C_+(0,1)$，限制其大于零，密度函数为
$\pi(\lambda)=\dfrac{2}{\pi}\cdot\dfrac{1}{1+\lambda^2},\quad\lambda>0.$
Local student’s t
把岭回归中的先验推广到每个自变量不同的场合，这样允许更多的波动性，即
$\begin{array}{rl}
\beta_j\mid\tau_j^2,\sigma^2\sim&\!\!\!\!N(0,\sigma^2\tau_j^2),\quad j=1,\dots,p,\\
\tau_j^2\mid\nu,\lambda\sim&\!\!\!\!\Gamma^{-1}\left(\dfrac{\nu}{2},\dfrac{\nu}{2\lambda}\right),\quad j=1,\dots,p,\\
\lambda\sim&\!\!\!\!C_+(0,1).
\end{array}$
这实际上就是 $NIG\left(0,\sigma^2,\dfrac{\nu}{2},\dfrac{\nu}{2\lambda}\right)$，因此边际为 $t$ 分布
$\beta_j\mid\sigma^2,\lambda\sim t\left(\nu,0,\dfrac{\sigma^2}{\lambda}\right),\quad j=1,\dots,p.$
Lasso
贝叶斯Lasso是正态-指数混合系统，即
$\begin{array}{rl}
\beta_j\mid\tau_j^2,\sigma^2\sim&\!\!\!\!N(0,\sigma^2\tau_j^2),\quad j=1,\dots,p,\\
\tau_j^2\mid\lambda\sim&\!\!\!\!\text{Exp}\left(\dfrac{\lambda^2}{2}\right),\quad j=1,\dots,p,\\
\lambda\sim&\!\!\!\!C_+(0,1).
\end{array}$
此时将 $\tau_j^2$ 积分掉可得边际为Laplace分布
$\beta_j\mid\sigma^2,\lambda\sim\text{Laplace}\left(0,\dfrac{\sigma}{\lambda}\right),\quad j=1,\dots,p.$
Elastic net
弹性网络是同时使用1-范数和2-范数进行惩罚的结果，按贝叶斯的观点，它可以通过正态混合得到，即
$\begin{array}{rl}
\beta_j\mid\tau_j,\sigma^2,\lambda_2\sim&\!\!\!\!N\left(0,\dfrac{\tau_j-1}{\tau_j}\cdot\dfrac{\sigma^2}{\lambda_2}\right),\quad j=1,\dots,p,\\
\tau_j\mid\sigma^2,\lambda_1,\lambda_2\sim&\!\!\!\!\Gamma_{1+}\left(\dfrac{1}{2},\dfrac{8\lambda_2\sigma^2}{\lambda_1^2}\right),\quad j=1,\dots,p,\\
\lambda_1\sim&\!\!\!\!C_+(0,1),\\
\lambda_2\sim&\!\!\!\!C_+(0,1),
\end{array}$
其中 $\Gamma_{1+}$ 表示支撑范围限制在 $(1,+\infty)$ 的截断伽马分布。将 $\tau_j$ 积分掉就可以得到表格中的条件先验密度形式。
Group Lasso
Group Lasso也可以通过正态混合得到，即
$\begin{array}{rl}
\pmb\beta_g\mid\tau_g^2,\sigma^2\sim&\!\!\!\!N\left(\pmb 0,\sigma^2\tau_g^2\pmb I_{m_g}\right),\quad g=1,\dots,G,\\
\tau_g^2\mid\lambda\sim&\!\!\!\!\Gamma\left(\dfrac{m_g+1}{2},\dfrac{\lambda^2}{2}\right),\quad g=1,\dots,G,\\
\lambda\sim&\!\!\!\!C_+(0,1).
\end{array}$
将 $\tau_g^2$ 积分掉就可以得到表格中的条件先验密度形式。
HyperLasso
HyperLasso是一种Adaptive Lasso的贝叶斯方法，其通过如下正态混合得到：
$\begin{array}{rl}
\beta_j\mid\phi_j^2\sim&\!\!\!\!N(0,\phi_j^2),\quad j=1,\dots,p,\\
\phi_j^2\mid\tau_j\sim&\!\!\!\!\text{Exp}(\tau_j),\quad j=1,\dots,p,\\
\tau_j\mid\nu,\lambda\sim&\!\!\!\!\Gamma\left(\nu,\dfrac{1}{\lambda^2}\right),\quad j=1,\dots,p,\\
\lambda\sim&\!\!\!\!C_+(0,1).
\end{array}$
按前面的正态-指数混合等价于Laplace先验知，上述混合等价于
$\begin{array}{rl}
\beta_j\mid\tau_j\sim&\!\!\!\!\text{Laplace}\left(0,\sqrt{2\tau_j}\right),\quad j=1,\dots,p,\\
\tau_j\mid\nu,\lambda\sim&\!\!\!\!\Gamma\left(\nu,\dfrac{1}{\lambda^2}\right),\quad j=1,\dots,p,\\
\lambda\sim&\!\!\!\!C_+(0,1).
\end{array}$
将 $\tau_j$ 积分掉就可以得到表格中的条件先验密度形式。
Horseshoe
Horseshoe先验是一种流行的贝叶斯惩罚先验，其混合模式为
$\begin{array}{rl}
\beta_j\mid\tau_j\sim&\!\!\!\!N(0,\tau_j^2),\quad j=1,\dots,p,\\
\tau_j\mid\lambda\sim&\!\!\!\!C_+(0,\lambda),\quad j=1,\dots,p,\\
\lambda\mid\sigma\sim&\!\!\!\!C_+(0,\sigma).
\end{array}$
由于半柯西先验可以表示为逆伽马和伽马密度混合，因此上述Horseshoe先验可以等价表示为：
$\begin{array}{rl}
\beta_j\mid\tau_j^2\sim&\!\!\!\!N(0,\tau_j^2),\quad j=1,\dots,p,\\
\tau_j^2\mid\omega\sim&\!\!\!\!\Gamma^{-1}\left(\dfrac{1}{2},\omega\right),\quad j=1,\dots,p,\\
\omega\mid\lambda^2\sim&\!\!\!\!\Gamma\left(\dfrac{1}{2},\lambda^2\right),\\
\lambda^2\mid\gamma\sim&\!\!\!\!\Gamma^{-1}\left(\dfrac{1}{2},\gamma\right),\\
\gamma\mid\sigma^2\sim&\!\!\!\!\Gamma\left(\dfrac{1}{2},\sigma^2\right).\\
\end{array}$
Horseshoe先验对大的系数向零压缩很小，尽管这一性质在理论上很好，但在实际应用中可能会出问题，特别是当参数是弱可识别的情形。此时，回归系数的后验均值估计可能并不存在，而且即使存在，Horseshoe先验也会产生不稳定低的MCMC抽样。
为此，Piironen and Vehtari（2017） 提出Regularized Horseshoe先验：
$\begin{array}{rl}
\beta_j\mid\tilde{\tau}_j^2,\lambda\sim&\!\!\!\!N(0,\tilde{\tau}_j^2\lambda),\quad j=1,\dots,p,\\
\lambda\mid\lambda_0^2\sim&\!\!\!\!C_+(0,\lambda_0^2),\\
\tau_j\sim&\!\!\!\!C_+(0,1),\quad j=1,\dots,p,\\
c^2\mid\nu,s^2\sim&\!\!\!\!\Gamma^{-1}\left(\dfrac{\nu}{2},\dfrac{\nu s^2}{2}\right).
\end{array}$
其中 $\tilde{\tau}_j^2=\dfrac{c^2\tau_j^2}{c^2+\lambda^2\tau_j^2}$，$\lambda_0=\dfrac{p_0}{p-p_0}\dfrac{\sigma}{\sqrt{n}}$，$p_0$ 表示对实际变量个数的先验猜测。
Spike-and-slab
这种方法使用一个质量集中在0附近的先验（the spike）与一个模糊先验（the slab）进行混合：
$\begin{array}{rl}
\beta_j\mid\gamma_j,\tau_j^2,\phi_j^2\sim&\!\!\!\!\gamma_j N(0,\tau_j^2)+(1-\gamma_j)N(0,\phi_j^2),\quad j=1,\dots,p,\\
\tau_j^2\sim&\!\!\!\!\Gamma^{-1}(0.5,0.5),\quad j=1,\dots,p,
\end{array}$
其中通过对 $\tau_j^2$ 取一个模糊先验使得可以基于数据估计the slab的方差。$\phi_j^2$ 固定为一个比较小的正数（比如0.001）来创建the spike。通过对 $\tau_j^2$ 取逆伽马先验，the slab部分的边际分布为柯西分布。
总结
上述这些方法可以通过R包 bayesreg 来进行。
- 在 $p<n$ 时，不同的惩罚先验效果类似于经典的惩罚方法效果。尽管经典的惩罚方法运行比贝叶斯方法快，但是它们没有给出一个精确的不确定性度量，所得的Bootstrap置信区间一般比贝叶斯可信区间要小很多。
- 在 $p>n$ 时，（Regularized）Horseshoe和HyperLasso先验在PMSE准则下要比绝大数其他先验表现好很多。
后记
本文内容参考自中国科学技术大学《贝叶斯分析》（2022Fall）课件。
如对本文内容有任何疑问，请联系 watthu@mail.ustc.edu.cn。
