从pymc3中推断的参数生成预测

Question

我遇到了一个常见问题，想知道是否有人可以帮忙。我经常想以两种模式使用pymc3：训练（即实际运行参数推断）和评估（即使用推断的参数生成预测）。

通常情况下，我希望得到预测的后验概率分布，而不仅仅是点估计（这是贝叶斯框架的优势之一，不是吗？）。当您的训练数据固定时，通常通过添加与观察变量类似的模拟变量来实现此目的。例如：

from pymc3 import *

with basic_model:

    # Priors for unknown model parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
    Y_sim = Normal('Y_sim', mu=mu, sd=sigma, shape=len(X1))

    start = find_MAP()
    step = NUTS(scaling=start)
    trace = sample(2000, step, start=start)

但是如果我的数据变化了怎么办？比如说我想基于新数据生成预测，但不想再重新运行推断。理想情况下，我希望有一个函数像 predict_posterior(X1_new, X2_new, 'Y_sim', trace=trace) 或者 predict_point(X1_new, X2_new, 'Y_sim', vals=trace[-1])， 可以直接将新数据运行到Theano计算图中。

我想我的问题的一部分与pymc3如何实现Theano计算图有关。我注意到函数model.Y_sim.eval似乎类似于我想要的，但它需要Y_sim作为输入，似乎只返回你给它的东西。

我想这个过程非常普遍，但我似乎找不到任何方法来做到这一点。非常感谢任何帮助。（还要注意，我有一个在pymc2中进行此操作的hack；由于Theano，pymc3中更难）

- santon

---

你在谈论从后验预测分布中进行抽样，看起来你做得很正确。不过我不确定你所说的“基于新数据”具体指什么。你是在讨论使用此分析的后验作为先验，以便基于额外数据进行推断吗？ - Chris Fonnesbeck

@ChrisFonnesbeck 我也对此很感兴趣，因为我们得到的后验概率是以迹的形式呈现的，我们无法使用它们来指定示例语法中的先验概率。 - recluze

twiecki在pymc3 gitter页面上向我指出了这个页面，它似乎解决了我面临的问题。我需要花些时间去理解所做的事情，但看起来很有前途。 - santon

有足够的数据，后验概率往往是（多元）正态分布。一个简单的方法是提取后验的均值和标准差，并用它来参数化正态先验进行后续分析。 - Chris Fonnesbeck

Answer 1

注意：此功能现已作为“pymc.sample_ppc”方法纳入核心代码中。请查看文档获取更多信息。

基于由twiecki发送给我的链接（截至2017年7月已失效），有几个技巧可以解决我的问题。第一个是将训练数据放入共享的theano变量中。这样可以在不破坏theano计算图的情况下稍后更改数据。

X1_shared = theano.shared(X1)
X2_shared = theano.shared(X2)

下一步，按照惯例构建模型并运行推断，但使用共享变量。

with basic_model:

    # Priors for unknown model parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1_shared + beta[1]*X2_shared

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

    start = find_MAP()
    step = NUTS(scaling=start)
    trace = sample(2000, step, start=start)

最后，正在开发一项功能（可能最终会添加到pymc3中），可以预测新数据的后验分布。

from collections import defaultdict

def run_ppc(trace, samples=100, model=None):
    """Generate Posterior Predictive samples from a model given a trace.
    """
    if model is None:
         model = pm.modelcontext(model)

    ppc = defaultdict(list)
    for idx in np.random.randint(0, len(trace), samples):
        param = trace[idx]
        for obs in model.observed_RVs:
            ppc[obs.name].append(obs.distribution.random(point=param))

    return ppc

接下来，传入您想要运行预测的新数据：

X1_shared.set_value(X1_new)
X2_shared.set_value(X2_new)

最后，您可以为新数据生成后验预测样本。

ppc = run_ppc(trace, model=model, samples=200)

变量ppc是一个字典，其中包含模型中每个观测变量的键。因此，在这种情况下，ppc['Y_obs']将包含一个数组列表，其中每个数组都是使用来自跟踪的单个参数集生成的。

请注意，您甚至可以修改从跟踪中提取的参数。例如，我有一个使用GaussianRandomWalk变量的模型，并且我想将预测生成到未来。虽然您可以允许pymc3对未来进行采样（即允许随机行走变量发散），但我只想使用对应于最后推断值的系数的固定值。这种逻辑可以在run_ppc函数中实现。

值得一提的是，run_ppc函数非常慢。它需要与运行实际推理一样多的时间。我怀疑这与theano的使用方式有关的某些低效率有关。

编辑：原链接似乎已失效。

- santon

Answer 2

@santon给出的回答是正确的，我只是进行了补充。

现在你不需要再编写自己的方法run_ppc。 pymc3提供了sample_posterior_predictive方法，它可以执行相同的操作。

- Ashok Rayal