22

我遇到了一个常见问题，我想知道是否有人可以提供帮助。我经常想在两种模式下使用 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 中有一个技巧可以做到这一点；由于 theano，在 pymc3 中更难。）

4

2 回答 2

13

注意：此功能现在作为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']将包含一个数组列表，每个数组都是使用来自 trace 的一组参数生成的。

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

还值得一提的是，该run_ppc功能非常缓慢。它花费的时间与运行实际推理的时间差不多。我怀疑这与使用 theano 的一些低效率有关。

编辑：最初包含的链接似乎已失效。

于 2015-10-22T21:16:10.063 回答
5

@santon 的上述回答是正确的。我只是补充一点。

现在您不需要编写自己的方法run_ppc。pymc3提供sample_posterior_predictive了同样的方法。

于 2019-04-26T09:19:53.370 回答