【Note】 贝叶斯数据分析(三)——层次化模型

前面我们介绍了贝叶斯框架下的单参数模型和多参数模型，但有些时候我们的模型并不是平的，而是具有层次化结构的——参数的影响一层层传递下来最后影响到观测。本章的内容十分深刻：这使得基于贝叶斯网络的模型成为可能。

引入

这里选用教材里给出的一个例子。我们有了70组历史实验数据，每组有\(n_i\)只小鼠，其中\(y_i\)只得了肿瘤。现在我们有了第71组数据——14只小鼠中有4只得了肿瘤，我们想知道这组小鼠得肿瘤的概率。（当然，我们不是去统计整体的得肿瘤的概率，因为组与组之间因为培养环境、操作人员等因素，其得肿瘤的概率可能会有差别）

事实上这是一个相当简单的问题——我们很容易就将其建模为之前研究过的单参数的二项分布模型（参见这里）。我们假设得肿瘤的概率服从\(\mathrm{Beta}(\alpha,\beta)\)，因此我们的后验分布即为\(\mathrm{Beta}(\alpha+4,\beta+10)\)

但我们如何得到超参数\(\alpha,\beta\)呢？前面提到了我们可以使用样本矩求出超参数——也就是矩估计。具体说来，我们求得70个肿瘤概率的均值方差.

根据\(\mathrm{Beta}\)分布有均值和方差

\[\mathbb{E}(\theta)=\frac{\alpha}{\alpha+\beta},\,\mathrm{Var}(\theta)=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}\]

求解得

\[\alpha=\frac{\mathbb{E}(\theta)^2(1-\mathbb{E}(\theta))}{\mathrm{Var}(\theta)}-\mathbb{E}(\theta),\,\beta=\frac{1-\mathbb{E}(\theta)}{\mathbb{E}(\theta)}\alpha\]

由此我们可以得到参数的后验分布，也就可以根据其最值、中位数、众数等做推断。

计算的R代码如下，数据文件可以在数据 下载到。

1
2
3
4
5

rats <- read.table('rats_data.txt',header=T)
E <- mean(rats$y[-71]/rats$n[-71]) # 0.1360065
s <- sd(rats$y[-71]/rats$n[-71]) # 0.1034923
alpha <- E^2*(1-E)/s^2-E # 1.356149
beta <- (1-E)/E*alpha # 8.615058

回顾

让我们重新审视下我们的模型。我们假设肿瘤率\(\theta\)服从一个超先验\(\mathrm{Beta}(\alpha,\beta)\)，每组实验都会从这个超先验中采样出一个肿瘤率，然后以这个肿瘤率为参数，从二项分布中采样得到观测到的数据\(y,n\)

事实上，这样的模型可以视为多参数模型的一个特例。我们的先验函数为\(p(\alpha,\beta,\theta)=p(\alpha,\beta)p(\theta\lvert\alpha\beta)\)，后验则相应地为\(p(\alpha,\beta,\theta\lvert y)\propto p(\alpha,\beta,\theta)p(y\lvert\alpha,\beta,\theta)=p(\alpha,\beta,\theta)\boldsymbol{p(y\lvert\theta)}\)。加粗部分说明了其和原始多参数模型的区别——观测\(y\)不再和参数\(\alpha,\beta\)直接相关，而是经过了\(\theta\)作为一个桥梁。

那么引入层次化模型有什么好处呢？

• 将模型分解为更小的子模块，每个子模块都更容易被理解，正如上面将整个模型拆成了两个易于理解的简单模型。
• 非层次化模型中，过少的参数往往难以很好的拟合模型
• 非层次化模型中，过多的参数很容易过拟合模型（因为有太多的条件概率）
• 在层次化模型中，有时甚至参数的个数可以多于数据的个数，仍能取得较合理的效果

框架

因此，层次模型的基本框架可以视为

\[y\lvert\theta,\phi\sim p(y\lvert\theta)\] \[\theta\lvert\phi\sim p(\theta\lvert\phi)\] \[\phi\sim p(\phi)\]

则我们的联合先验和后验分别为

\[p(\theta,\phi)=p(\phi)p(\theta\lvert\phi)\] \[p(\theta,\phi\lvert y)\propto p(\phi)p(\theta\lvert\phi)p(y\lvert\theta,\phi)=p(\phi)p(\theta\lvert\phi)p(y\lvert\theta)\]

通常我们感兴趣的量为\(p(\phi\lvert y)\)，这可以用两种方法求得，各有优劣

• 直接积分法\(\int p(\theta,\phi\lvert y)\,\mathrm{d}\theta\)，但积分操作在高维空间里并不容易求得。
• 条件密度法\(p(\theta,\phi\lvert y)/p(\theta\lvert\phi,y)\)，但分母的归一化因子通常是与\(\phi\)有关的，有可能其计算并不简单

在层次化模型中，超先验通常会选为无信息分布——因为它和我们的观测“较远”，一般我们对其没有太多的了解。然而，超先验的选取通常是十分重要的，主要原因是一个其有可能导致一个improper的后验，因此我们必须加以选择。接下来的例子就会阐释这一点。

回到最初的例子——二项分布

求解后验分布

在之前的例子中，我们是使用前70个样本给出了超参数的一个估计。现在我们给超参数赋予一个先验，然后将71个样本都输入，来看\(\alpha,\beta\)的后验分布。

联合后验为

\[p(\alpha,\beta,\theta\lvert y)\propto p(\alpha,\beta) \prod_{J=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta_j^{\alpha-1}(1-\theta_j)^{\beta-1}\prod_{J=1}^{J}\theta_j^{y_j}(1-\theta_j)^{n_j-y_j}\]

\(\theta\)的条件后验为一个\(\mathrm{Beta}\)分布，为

\[p(\theta\lvert \alpha,\beta,y)=\prod_{J=1}^{J}\frac{\Gamma(\alpha+\beta+n_j)}{\Gamma(\alpha+y_j)\Gamma(\beta+n_j-y_j)}\theta_j^{\alpha+y_j-1}(1-\theta_j)^{\beta-1+n_j-y_j}\]

此时超先验的后验分布（边缘）为

\[p(\alpha,\beta\lvert y)\propto p(\alpha,\beta)\prod_{J=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\frac{\Gamma(\alpha+y_j)\Gamma(\beta+n_j-y_j)}{\Gamma(\alpha+\beta+n_j)}\]

超先验的选择

为了后验

因此，我们选取的超先验\(p(\alpha,\beta)\)需要使得上式为一个proper的后验。但要注意到，上式的连乘项当\(\alpha,\beta\)趋于无穷时其值趋于1，因此其积分显然不收敛，因此\(p(\alpha,\beta)\)在无穷远处必须（比较迅速地）趋于零才有可能将这个积分“压为”有限值，或者说其必须为轻尾(light tail)的。

一个很常用的变换是将其变换为\((\log(\alpha/\beta),\log(\alpha+\beta))\)，这么做的一个原因是，\(\alpha+\beta\)表示着先验带来的“数据量”，而\(\log(\alpha+\beta)=\mathrm{logit}(\alpha/(\alpha+\beta))\)，而\(\alpha/(\alpha+\beta)\)是先验带来的“样本均值”。

然而，很遗憾地，可以证明其仍然不能使后验变得proper的(提示：考虑\(\alpha/\beta\)为非零常数并取极限\(\alpha+\beta\rightarrow\infty\))

基于上面的思想，再尝试不同的形式，可以发现\(p(\frac{\alpha}{\alpha+\beta},\alpha+\beta)\propto 1\)仍然是improper的，而\(p(\frac{\alpha}{\alpha+\beta},(\alpha+\beta)^{-1/2})\propto 1\)是可以采用的，此时等价的，为\(p(\alpha,\beta)\propto (\alpha+\beta)^{-5/2}\)，或写为\(p(\log(\alpha/\beta),\log(\alpha+\beta))\propto \alpha\beta(\alpha+\beta)^{-5/2}\)。

当然，读者也可以选择其他的先验，只要其使得后验proper就可以。

为了简单

事实上，一个更加偷懒的方式就是，我们不去选择一个跨越全空间的均匀分布，而只取其中一个测度有限的子集，这样一来我们的先验就不存在improper的问题，后验也自然是proper的。当然，一个不可避免的问题就是，如果子集选择得不合理，可能导致整个推断都将是失败的。然而借助于计算机的力量，我们可以比较轻松的避免这个问题。

仍然使用刚才的例子，我们通过矩估计其实得到了超参数\(\alpha,\beta\)的大致取值范围。因此我们不妨认为\(\alpha\)取值在\(0-5\)中间，\(\beta\)取值在\(0-20\)中间且为均匀分布。那么这样的后验分布会是什么样呢？

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

rats <- read.table('rats_data.txt',header=T)
n <- rats$n
y <- rats$y
GRID.SIZE <- 200

log.post.cut.ab <- function(A,B) {
    Vectorize(function(a,b) {
        if (a > A || b > B) return (-Inf)
        sum(lbeta(a+y,b+n-y) - lbeta(a,b))
    })
}

test <- function(A, B) {
    ALPHA <- (1:GRID.SIZE) / (GRID.SIZE / A)
    BETA <- (1:GRID.SIZE) / (GRID.SIZE / B)
    post.ab.result <- outer(ALPHA, Y=BETA, FUN=log.post.cut.ab(A,B))
    post.ab.result <- exp(post.ab.result - max(post.ab.result))
    contour(ALPHA, BETA, post.ab.result)
}

test(5,20)

可以看到其后验分布被“切断了”，这提示我们在\(\beta\)上的选取不够合理。经过调整后，test(5,40)的效果相对较好。

与之对比，我们做出用我们上一小节求出的一个符合要求的先验，并同样对其后验分布作图。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

log.post.ab <- Vectorize(function(a,b) {
    -2.5 * log(a+b) + sum(lbeta(a+y,b+n-y) - lbeta(a,b))
})

test1 <- function() {
    ALPHA <- (1:GRID.SIZE) / (GRID.SIZE / 5)
    BETA <- (1:GRID.SIZE) / (GRID.SIZE / 40)
    post.ab.result <- outer(ALPHA,Y=BETA, FUN=log.post.ab)
    post.ab.result <- exp(post.ab.result-max(post.ab.result))
    contour(ALPHA, BETA, post.ab.result)
}

par(mfrow=c(1,2))
test(5,40)
test1()

可以看到二者的形状基本相似，其差别在于左边的中心相对来说更靠右上方一些，也就是\(\alpha,\beta\)均会略大。这是十分自然的，因为左边的先验是均匀分布，而右边的先验是\((\alpha+\beta)^{-5/2}\)，这个先验更加青睐小的\(\alpha+\beta\)。

十分重要的例子——正态分布

模型假设

一共有\(J\)组，第\(j\)组有\(n_j\)个样本。组内服从均值\(\theta_j\)未知，方差\(\sigma^2\)已知且共享的正态分布。如果读者对方差分析有所了解，对这个模型应该不陌生。方差分析我在Note-Analysis-of-Variance也有介绍。而这里我们假设各组正态分布的均值是独立同分布于一个正态分布的，即\(p(\theta_1,\ldots,\theta_J)=\prod_{j=1}^{J}\mathcal{N}(\theta_j\lvert\mu,\tau^2)\)

可以求出联合后验为

\[p(\theta,\mu,\tau\lvert y)\propto p(\mu,\tau)\prod_{j=1}^{J}\mathcal{N}(\theta_j\lvert\mu,\tau^2)\prod_{j=1}^{J}\mathcal{N}(\bar{y}_{\cdot j}\lvert\theta_j,\sigma_j^2)\]

其中\(\bar{y}_{\cdot j}=\sum_{i=1}^{n_j}y_{ij}/n_j,\,\sigma_j^2=\sigma^2/n_j\)。进一步求得\(\theta\)的条件后验分布以及超参数\(\mu,\tau\)的边缘后验分布为

\[p(\theta_j\lvert\mu,\tau,y)\sim\mathcal{N}(\widehat{\theta}_j,V_j)\qquad\widehat{\theta}_j=\frac{\frac{\bar{y}_{\cdot j}}{\sigma_j^2}+\frac{\mu}{\tau^2}}{\frac{1}{\sigma_j^2}+\frac{1}{\tau^2}},\,V_j=\frac{1}{\frac{1}{\sigma_j^2}+\frac{1}{\tau^2}}\] \[p(\mu,\tau\lvert y)\propto p(\mu,\tau)\prod_{j=1}^{J}\mathcal{N}(\bar{y}_{\cdot j}\lvert\mu,\sigma_j^2+\tau^2)\]

在这里，我们选用无信息先验\(p(\mu,\tau)\propto p(\mu\lvert\tau)p(\tau)\propto p(\tau)\)，再尝试决定\(p(\tau)的形式使得\)p(y)$合理。

注意到\(p(\mu,\tau\lvert y)=p(\mu\lvert\tau,y)p(\tau\lvert y)\)，我们首先可以算出

\[p(\mu\lvert\tau,y)\sim\mathcal{N}(\widehat{\mu},V_\mu),\quad \widehat{\mu}=\frac{\sum_{j=1}^{J}\frac{\bar{y}_{\cdot j}}{\sigma_j^2+\tau^2}}{\sum_{j=1}^{J}\frac{1}{\sigma_j^2+\tau^2}},\,V_\mu^{-1}=\sum_{j=1}^{J}\frac{1}{\sigma_j^2+\tau^2}\]

由此可得

\[p(\tau\lvert y)\propto\frac{p(\mu,\tau\lvert y)}{p(\mu\lvert\tau,y)}\propto \frac{p(\tau)\prod_{j=1}^J\mathcal{N}(\widehat{y}_{\cdot j}\lvert\mu,\sigma_j^2+\tau^2)}{\mathcal{N}(\mu\lvert\widehat{\mu},V_\mu)}\]

注意到上式对于任意\(\mu\)都是成立的，因此我们不妨取\(\mu=\widehat{\mu}\)简化计算，得

\[p(\tau\lvert y)\propto p(\tau)V_\mu^{1/2}\prod_{j=1}^{J}(\sigma_j^2+\tau^2)^{-1/2}\exp\left\{-\frac{(\bar{y}_{\cdot j}-\widehat{\mu})^2}{2(\sigma_j^2+\tau^2)}\right\}\]

可以证明，\(p(\tau)\propto 1\)是合适的，逆卡方分布也是合适的，而\(p(\log\tau)\propto 1\)是不合适的

实际计算

因此实际使用中我们使用抽样的方法分析其联合后验。抽样步骤如下

1. 从\(p(\tau\lvert y)\)中抽样出\(\tau_k\)
2. 从\(p(\mu\lvert \tau_k,y)\)中抽样出\(\mu_k\)
3. 从\(p(\theta\lvert \mu_k,\tau_k,y)\)中抽样出\(\theta_k\)（事实上，由于\(\theta_1,\ldots,\theta_J\)彼此条件独立，我们这一步可以直接抽样出\(J\)个样本\(\theta_{1,k},\ldots,\theta_{J,k}\)。

如果不用贝叶斯方法？

我们之前讲过，还有一种确定先验的方式是利用样本信息做点估计。然而在这个例子中，求取\(\mu,\tau^2\)的无偏估计得到

\[\widehat{\mu}=\bar{y}_{\cdot\cdot},\,\widehat{\tau}^2=(\mathrm{MS_B}-\mathrm{MS_W})/n\]

其中\(\mathrm{MS_B},\mathrm{MS_W}\)分别表示组内和组间的残差平均值。但这样做并不好，一个显而易见的问题就是\(\widehat{\tau}^2\)有可能为负值。同时这样的求法意味着我们丢弃了一些\(\theta_j\)的不确定性。

例子分析

现在有八家培训机构，我们想知道这些培训机构是否真的能提高被试者的成绩，因此我们进行了一系列实验，得到了成绩提升的均值和方差。我们想知道\(A\)的异军突起是否科学？还是只是一次偶然？

我们先考察两种不使用贝叶斯模型的方法来求解这一问题

独立预测

也就是认为这八组正态分布彼此独立没有联系。那么\(A\)的效果就是均值为\(28\)，标准差为\(14.9\)的正态分布。

在此假设下，A的效果比\(28\)更好的概率是\(50\%\)。

汇集预测

这里我们使用\(\bar{y}_{\cdot\cdot}\)估计\(\theta\)，求出后验均值和后验方差

\[\bar{y}_{\cdot\cdot}=\frac{\sum_{j=1}^{J}\bar{y}_{\cdot j}/\sigma_j^2}{\sum_{j=1}^{J}1/\sigma_j^2}=7.69\] \[(\sum_{j=1}^{J}1/\sigma_j^2)^{-1/2}=4.07\]

1
2
3
4
5

y<-c(28,8,-3,7,-1,1,18,12)
s<-c(15,10,16,11,9,11,10,18)

pool.mean <- sum(y/s^2) / sum(1/s^2)
pool.s <- sqrt(1 / sum(1/s^2))

那么\(A\)的效果就是均值为\(7.7\)，标准差为\(4.1\)的正态分布。构造出其\(95\%\)后验置信区间约为\([8\pm8]\)

在此假设下，A的效果比\(7.7\)更差的概率是\(50\%\)。且A的效果比最差的C更差的概率也是\(50\%\)

可以看到，上面这两个模型都不能给出满意的结果。我们更想把两者的信息结合起来。一个可行的做法是认为每组的均值为\(\widehat{\theta}_j=\lambda_j\bar{y}_{\cdot j}+(1-\lambda_j)\bar{y}_{\cdot\cdot}\)，也就是做一个加权平均，但\(\lambda_j\)的选择也需要仔细考量。

贝叶斯模型

我们不妨取\(p(\tau)\propto 1\)来求取后验分布，绘图如下

1
2
3
4
5
6

post.tau.given.y <- Vectorize(function(tau) {
  mu.hat <- sum(y/(s^2+tau^2)) / sum(1/(s^2+tau^2))
  mu.V <- 1 / sum(1/(s^2+tau^2))
  mu.V^(0.5)*prod((s^2+tau^2)^(-0.5)*exp(-(y-mu.hat)^2/(2*(s^2+tau^2))))
})
curve(post.tau.given.y(x), 0, 30)

可以看到后验会倾向于较小的\(\tau\)。

接着我们考察条件后验分布\(\mathbb{E}(\theta_j\lvert\tau,y)\)

1
2
3
4
5
6
7
8
9
10
11

post.thetamean.given.tau.y <- function(id) {
  Vectorize(function(tau) {
    mu.hat <- sum(y/(s^2+tau^2)) / sum(1/(s^2+tau^2))
    (y[id]/s[id]^2+mu.hat/tau^2) / (1/s[id]^2+1/tau^2)
  })
}
plot(0, type='n', xlim=c(0,30), ylim=c(-5,25))
for (i in 1:8) {
  curve(post.thetamean.given.tau.y(i)(x), 0, 30, add=T)
  text(20+rnorm(1), post.thetamean.given.tau.y(i)(20)+0.5, LETTERS[i])
}

而关于\(A\)的效果的回答，我们可以通过上述的抽样流程，从后验分布中抽取\(\theta_1\)的，然后通过其分布的信息（如分位数）就可以给出其大于\(28\)的概率(p-value)。甚至只需简单地统计10000个样本中有多少个超过\(28\)即可。但关于如何从这样的一个分布中抽样，这里咱不提及，在后面关于贝叶斯计算的章节将会详细地讨论这一问题。