从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中更难.)

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']将包含一个数组列表,每个数组都是使用trace中的一组参数生成的.

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

值得一提的是,该run_ppc功能非常慢.它需要与运行实际推理一样多的时间.我怀疑这与使用theano的一些低效率有关.

编辑:最初包含的链接似乎已经死了.

Answer 2

@santon 的上述答案是正确的。我只是补充一下。

现在您不需要编写自己的方法run_ppc。pymc3提供了sample_posterior_predictive执行相同操作的方法。