linear-regression - 使用 PyMC3 进行基本贝叶斯线性回归预测

Question

我想使用我的 PyMC3 LR 模型在新数据可用时获得预测变量 y 值的 80% HPD 范围。因此，为 x 的新值推断出 y 值的可靠分布，而不是在我的原始数据集中。

型号:

with pm.Model() as model_tlr:
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10)
    epsilon = pm.Uniform('epsilon', 0, 25)

    nu = pm.Deterministic('nu', pm.Exponential('nu_', 1/29) + 1)
    mu = pm.Deterministic('mu', alpha + beta * x)

    yl = pm.StudentT('yl', mu=mu, sd=epsilon, nu=nu, observed=y)

    trace_tlr = pm.sample(50000, njobs=3)

在 burnin 之后，我从后部采样并获得 HPD

ppc_tlr = pm.sample_ppc(btrace_tlr, samples=10000, model=model_tlr)
ys = ppc_tlr['yl']
y_hpd = pm.stats.hpd(ys, alpha=0.2)

这对于可视化集中趋势周围的 HPD 非常有用(使用 fill_between)

但我现在想使用该模型在 x=126.2 时获取 y 的 HPD(例如)并且初始数据集不包含观察到的x=126.2

我对后验采样的理解是，数据集中每个可用的 x 值都有 10k 个样本，因此 ys x=126.2 的，因为它没有被观察到。

基本上，有没有一种方法可以使用我的模型从预测值 x=126.2 中获取可信值的分布(基于模型)，该预测值仅在模型构建后才可用？如果是，怎么办？

谢谢

编辑:
找到SO Post其中提到

Function under development (will likely eventually get added to pymc3) that will allow to predict posteriors for new data.

这个存在吗？

Answer 1

好的，所以这是可能的，或多或少如上述 SO 帖子中所述。但是，此后 PyMC3 中添加了一个 sample_ppc 函数，这使得作者的 run_ppc 变得多余。

首先，为 x 设置一个 Theano 共享变量。

from theano import shared
x_shared = shared(x)

然后在构建模型时使用 x_shared。

模型建立后，添加新数据并更新共享变量

x_updated = np.append(x, 126.2)
x_shared.set_value(x_updated)

使用原始跟踪和模型对象重新运行 PPC 样本生成器

new_ppc = pm.sample_ppc(btrace_tlr, samples=10000, model=model_tlr)

找到新数据的后验样本

sample = new_ppc['yl'][:,-1]

然后我可以得到 HPD

pm.stats.hpd(sample)

array([ 124.56126638, 128.63795388])

Sklearn 让我认为应该有一个简单的 predict 接口(interface)...