pymc 3.0 预测后验分布

Question

我正在将一个非常简单的示例从 pymc 2.3 转换为 pymc 3.0，并且似乎无法弄清楚如何从预测后验分布中采样（或获取 MAP）。按照文档中的建议（7.3 模型检查和诊断：拟合优度），我可以使用 pymc 2.3 通过添加未观察到的随机变量从该分布中采样。这是link to the notebook。一切似乎都运行良好。

但是，当我尝试在 pymc 3.0 中执行此操作时，发生了两件奇怪的事情。

1. MAP 值甚至不接近，好像未观察到的随机性影响了最小化，未观察到的随机性的 MAP 值是错误的？
2. NUTS 采样器不会改变未观察到的随机变量的值，因此轨迹只是一个单一的值 10。

当然，如果观察到的随机指标不存在，这个值是有意义的。在 pymc 3.0 中，为预测后验查找 MAP 值和采样有何变化？

更新：

这是一个最小的例子，说明我做错了什么：

import pymc as mc
with mc.Model() as model:
    p = mc.Beta('p',2,2)
    surv_sim = mc.Binomial('surv_sim',n=20,p=p)
    surv = mc.Binomial('surv',n=20,p=p,observed=15)

with model:
    step = mc.step_methods.HamiltonianMC(vars=model.vars) #Again must specificy
                                                          #model.vars or else
                                                          #only continuous values will
                                                          #be sampled

with model:
    trace = mc.sample(100000,step)

mc.traceplot(trace);

with model:
    map_est = mc.find_MAP(vars=model.vars) #Must Specify model.vars, or else
                                           #only continuous stochastics will
                                           #be fit, however this fit will be 
                                           #horrible

我想我自己已经部分回答了这个问题。首先是抽样。据我所知，NUTS 采样器无法处理离散变量。然而，Metropolis 采样器肯定可以处理离散变量（已在 issue #235 中解决），HamiltonianMC 采样器似乎也可以处理离散变量。

然而，离散随机指标的 MAP 估计注定是糟糕的。即使您指定使用所有变量（不是默认值）。因为离散随机对数概率函数返回任何离散随机的底的对数概率，所以无论使用哪个 scipy 最小化函数，都会陷入局部最优。最小化只会优化连续随机变量，同时保持离散变量固定，因为离散变量中的任何小步骤都不会导致对数概率的提高。我不确定这是否是一个错误，或者仅仅是当您同时具有离散变量和连续变量时查找 MAP 估计值的基本限制。

更新 2：

如果您使用多个步骤，采样效果会更好。对于上述模型，您可以使用。

with model:
    step1 = mc.step_methods.NUTS(vars=[p])
    step2 = mc.step_methods.Metropolis(vars=[surv_sim])

with model:
    trace = mc.sample(100000,[step1,step2])

Answer 1

最近有很多关于 PyMC3 的开发活动，所以我想我会为你的问题发布一个小答案，尽管你一年多前自己就知道了。

首先，pymc3.find_MAP 函数现在尝试处理离散变量，因此要找到 MAP，您只需这样做：

import pymc3 as pm

m = pm.Model()
with m:
    p        = pm.Beta('p',2,2)
    surv_sim = pm.Binomial('surv_sim', n=20, p=p)
    surv     = pm.Binomial('surv', n=20, p=p, observed=15)

    map_est = pm.find_MAP()

map_est

对我来说，这会导致以下结果

{'p': array(0.6190476210094416), 'surv_sim': 10.0} 

这似乎不正确。也许这值得一个错误报告。

从后验分布采样现在需要明确选择步进方法，对于离散变量pymc3.step_methods.Metropolis 是可接受的选择。它也适用于本例中的连续变量：

with m:
    step1 = pm.step_methods.Metropolis(vars=[p])
    step2 = pm.step_methods.Metropolis(vars=[surv_sim])
    trace = pm.sample(10000, [step1, step2])