贝叶斯层次型模型参数估计 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)

超参数 α 即shape 为5 β 即scale为1 生成了5000个样本

先说说第一个问题

什么是贝叶斯参数估计

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

p(θ|y)=p(y|θ)p(θ)/p(y)∝p(y|θ)p(θ)

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

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

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

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

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

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

层次型模型(Hierarchical model)

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

p(y,θ,α,β)∝θye−θy!×βαΓ(α)θα−1e−βθ 这时候就可以把 α 和 β 的等比例概率化简出来了 p(α|y,θ,β)∝βαθαΓ(α) p(β|y,θ,α)∝βαe−βθ 这里我们默认 α 和 β 的先验是uniform的，所以直接省略了，当然我们也可以给 α 和 β 一个先验分布。 这样三个参数的条件概率就都写出来了，这时候就可以采样了。

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

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