贝叶斯统计/PyMC 学习笔记（一）

最近在读 Doing Bayesian Data Analysis，由于这本书的代码是用 R 写的，而我暂时也没有学 R 的打算，所以一边看书的同时一边学习用 PyMC3 来实现书中的例子/习题。这篇先写写 pymc 的使用。

一些有用的链接

PyMC3 repo
PyMC3 repo 提供的 example
Probabilistic Programming and Bayesian Methods for Hackers （一本挺好的开源书，这本书是用 pymc2，这里用了它提供的 matplotlib 样式）
Doing Bayesian Data Analysis - python/pymc3 version

基本用法

例如我们要测试一枚有偏的硬币，出现 head 的概率为 0.7。先把这个硬币抛上 100 遍:

flips = stats.bernoulli.rvs(p=0.7, size=100)
plt.hist(flips, bins=2)

假设我们对这个硬币一无所知，则先验概率为均匀分布。定义模型：

with pm.Model() as model:
    p = pm.Uniform('p', 0.1, 0.9)
    obv = pm.Bernoulli('flips', p, observed=flips) # 观察到的数据

然后定义取样的算法和初始点，就可以进行 sampling：

with model:
    start = pm.find_MAP(vars=[p])
    step = pm.NUTS() # 书中详细介绍了 Metropolis 方法，pymc 也是有提供的

    trace = pm.sample(2000, start=start, step=step, progressbar=False)

p 的后验分布：

burnin = 1000
thin = 10

pm.traceplot(trace[burnin::thin], ['p'])

pm.summary(trace[burnin::thin], ['p'])

p:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.670            0.051            0.005            [0.584, 0.765]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.580          0.639          0.673          0.703          0.762

然后还有自相关检验：

pm.autocorrplot(trace[burnin::thin], ['p'])

以上。