贝叶斯层次型模型参数估计 Bayesian hierarchical model parameter estimation with Stan

再来总结下贝叶斯参数估计，分为以下几部分：
1. 先说说贝叶斯参数估计
2. 再说说层次型模型，指的就是超参数（Hyper parameter）的选择
3. 用R+stan的Hamiltonian MC把这些参数（数据分布的参数和超参数）都采出来

这里我们用一个例子来演示怎么估计参数。
我们使用一个人工的数据，每天超市里一件商品的销售量y。
生成数据的R代码：

nSamples <- 5000;

theta <- rgamma(nSamples,5,1)

y <- rpois(nSamples, theta)

超参数$\alpha$即shape 为5 $\beta$即scale为1 生成了5000个样本

先说说第一个问题

什么是贝叶斯参数估计

先说说贝叶斯参数估计的时候用的贝叶斯定理：

$$p(\theta|y) = p(y|\theta)p(\theta)/p(y) \propto p(y|\theta)p(\theta)$$

这里$p(\theta|y)$是后验概率 $p(y|\theta)$是统计模型也叫似然 $p(\theta)$是先验概率。

参数估计估计的就是这个$\theta$，要区别的是如果求$p(y|\theta)$的最大值就是MLE(Maxium Likelihood Estimation)，求$p(\theta|y)$的最大值就是MAP (Max A Posterior)。MAP不能算是贝叶斯估计，因为它没有整合所有的不确定性。

本文的例子是一个计数的例子，应该使用泊松分布来作为似然。即：

$$y \sim poisson(\theta|n)$$
其中 $\theta$是泊松分布的参数，n是实验次数。然后是先验，这里选择gamma分布作为参数的先验，即：
$$\theta \sim gamma(\alpha,\beta)$$
其中 $\alpha$和 $\beta$就是超参数。选择gamma分布是以为poisson分布和gamma分布的共轭性。简单来说共轭就是说后验的分布跟先验的分布是一样的，只是参数有区别，而这个区别来自于似然。指数分布族的各种分布都有共轭分布，这给计算带来了很大的方便。

有了先验和似然，就可以推导后验分布了：

$$p(\theta|y) \propto p(y|\theta)p(\theta) = poisson(y|\theta)*gamma(\theta|\alpha,\beta)$$
$$＝ \frac{\theta^y e^{-\theta}}{y!} \times \frac{\beta^\alpha}{\Gamma(\alpha)}\theta^{\alpha-1}e^{-\beta\theta}$$
有了这个联合概率分布之后就可以去掉不相关的变量了，因为这里求的是 $\theta$, 所有跟它不相关的都可以去掉，不过等号不再成立，此处应该是等比例于。
$$p(\theta|y) \propto \theta^{(y+\alpha)-1}e^{-(1+\beta)\theta}$$
这个式子很像gamma分布的kernel。什么是kernel？gamma的密度函数可以分为两部分： $\frac{\beta^\alpha}{\Gamma(\alpha)}$和 $\theta^{\alpha-1}e^{-\beta\theta}$前面那一块叫normalize constant，后面的叫 kernel。也就是说后面这一块才是这个密度函数里面最管用的。
这样可以就可以把后验写成这个样了。
$$p(\theta|y) \propto gamma(y+\alpha, 1+\beta)$$
这样假设 $\alpha$和 $\beta$都知道的情况下，我们就可以根据y来采样 $\theta$的值了。这里可能有时候想不太明白为什么这个式子里面随便取一个值就是 $\theta$的值。上面这个式子里面右边 $p(\theta|y)$的意义就是 $\theta$在给定y时的概率。也正是我们要求的概率，可以把整个 $p(\theta|y)$看作一个整体，有助于理解
下面要求 $\alpha$和 $\beta$，就得了解下什么是层次型模型。

层次型模型(Hierarchical model)

层次型模型很简单就是说在数据这一层上面还有一层别的分布。这里就是指的这些超参$\alpha$和$\beta$。层次型模型不止两层，也可以是N层，但是层数太多了就会让整个模型效率变低。
怎么求超参数？先把刚才的联合概率密度函数写出来：

$$p(y,\theta,\alpha,\beta) \propto \frac{\theta^y e^{-\theta}}{y!} \times \frac{\beta^\alpha}{\Gamma(\alpha)}\theta^{\alpha-1}e^{-\beta\theta}$$
这时候就可以把 $\alpha$和 $\beta$的等比例概率化简出来了
$$p(\alpha|y,\theta,\beta) \propto \frac{\beta^\alpha\theta^\alpha}{\Gamma(\alpha)}$$
$$p(\beta|y,\theta,\alpha) \propto \beta^\alpha e^{-\beta\theta}$$
这里我们默认 $\alpha$和 $\beta$的先验是uniform的，所以直接省略了，当然我们也可以给 $\alpha$和 $\beta$一个先验分布。
这样三个参数的条件概率就都写出来了，这时候就可以采样了。

R+stan参数估计

stan是一个专门针对于概率的编程语言，非常智能，前面的公式推导可以基本上都省区，stan会自动计算。但是为了方便理解，还是写出来了。stan的代码如下：

data{
  int<lower=0> N;
  int<lower=0> y[N];
}
parameters{
  real theta;
  real<lower=0,upper=10> a;
  real<lower=0,upper=10> b;
}
model{
  theta ~ gamma(a,b);
  y ~ poisson(theta);
}

data块说明了用在这个model里面的数据，y就是刚才生成的数据。
parameters块里面有要估计的三个参数。
model块里面说明数据是怎么生成的，可以看出来这个是个层次型模型。
下面是执行的R语言

library(rstan)

dat <- list(N = 5000,
                    y = y)

fit <- stan(file = 'poisson.stan', data = dat,
            iter = 4000)

print(fit)

下面是输出的结果：

          mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
theta     5.05    0.00 0.03     4.99     5.03     5.05     5.07     5.11  4347    1
a         6.70    0.04 2.29     1.81     5.07     7.12     8.60     9.85  3331    1
b         1.51    0.01 0.69     0.35     1.01     1.46     1.94     3.06  3888    1
lp__  15609.88    0.02 1.29 15606.61 15609.35 15610.21 15610.80 15611.31  2826    1

从这个结果里面可以看到$\theta$的估计值是5，这个说明每天这个商品平均被卖出去的量是在5个左右，这个跟gamma(5,1)的期望是一样的。$\alpha$和$\beta$的估计值跟真实的参数有些偏差，但是因为我们本来就没有找两个参数的准确值，只是他们的等比例的值，此时这两个参数的比例已经跟真实参数非常相近了。