贝叶斯模型

要建立一个贝叶斯模型，我们需要确定模型的输出（outcome）符合什么分布，这意味着我们需要选择一个似然函数。

eg.

$$ height \sim Normal(\mu, \sigma) $$

如何选择似然函数

选择似然函数的时候，要考虑：

  • 数据本身的情况，比如是整数还是浮点数，是否有上下界等（Epistemological, which means philosophy behind knowledge）。
  • 要描述的对象在现实中是如何分布的，比如大量随机影响下，数据常常（天然地）呈正态分布（Ontological, which means the nature of the world）。

所选择的似然函数常常会需要传入参数(parameters)，对于每个参数，我们需要在一开始估计它的分布，这个分布称为这个参数的先验(prior)。

eg.

$$ \sigma \sim Exp(?) $$

有些参数可能和你所知道的其他输入有关，此时，这些输入称为自变量（predictor variable）。

eg.

$$ \mu = \alpha + \beta \times (weight - 150) $$

这里的 $\alpha$ 和 $\beta$ 是新增加的参数，而 $weight$ 就是自变量。

所以我们整个的模型就是：

$$ \displaylines { height & \sim & Normal(\mu, \sigma) \\ \mu & = & \alpha + \beta \times weight \\ \alpha & \sim & Normal(?, ?) \\ \beta & \sim & Normal(?, ?) \\ \sigma & \sim & Exp(?) } $$

现在我们需要确定模型中这些 $?$ 的值，并验证这些值是相对好的，这个过程称为 prior predictive check。

Prior predictive check

Prior predictive check 的目的是检查模型中的先验分布是否合理。

一般来说 prior predictive check 的方法是选定一组 prior，然后用这组 prior 生成一些数据，看这些数据是否合理。

比如以上面的例子来说，如果取以下的 prior：

$$ \displaylines { height & \sim & Normal(\mu, \sigma) \\ \mu & = & \alpha + \beta \times (weight - 150) \\ \alpha & \sim & Normal(0, 1) \\ \beta & \sim & Normal(0, 1) \\ \sigma & \sim & Exp(1) } $$

那么模型可以生成（大量的）负数的 height，这显然是不合理的。

根据经验，人的身高大约在 165 厘米左右，且身高和体重呈正相关，因此一组更好的 prior 可以是：

$$ \displaylines { height & \sim & Normal(\mu, \sigma) \\ \mu & = & \alpha + \beta \times (weight - 60) \\ \alpha & \sim & Normal(165, 20) \\ \beta & \sim & Normal(0.5, 0.2) \\ \sigma & \sim & Exp(1) } $$

当然你也可以搜集更多的知识选择更好的 prior，或者换用更合理的分布函数，比如让 $\beta \sim \text{log-Normal}(?, ?)$。

如何学习参数的值

一个常用的方法是使用 MCMC（Markov chain Monte Carlo）。

例如，我们对于以下模型：

$$ \displaylines { W & \sim & Binomial(9, p) \\ p & \sim & Normal(0.5, 0.1) } $$

可以使用 MCMC 方法来学习 $p$ 的值，使用 rethinking 包如下：

ulam(
    alist(
        W ~ dbinom(W + L, p), # binomial likelihood
        p ~ dnorm(0.5, 0.1) # uniform prior
    ),
    data = list(W = 6, L = 3)
)

功能上相当于：

W <- 6
L <- 3
# MCMC 的一个 chain
n_samples <- 1000
p <- rep(NA, n_samples)
# prior
p[1] <- rnorm(1, 0.5, 0.1)
for (i in 2:n_samples) {
    # 选取一个新的 p
    p_new <- rnorm(1, p[i - 1], 0.1)
    # 比较新旧 p 的结果
    q0 <- dbinom(W, W + L, p[i - 1])
    q1 <- dbinom(W, W + L, p_new)
    # 判断是否接受新的 p
    p[i] <- ifelse(runif(1) < q1 / q0, p_new, p[i - 1])
}