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

Question

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

型号：

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 post 所述。但是，此后 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)

新数据的post先验采样是

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

然后我可以通过

获得HPD

pm.stats.hpd(sample)

array([ 124.56126638, 128.63795388])

Sklearn 让我觉得应该有一个简单的 predict 界面...

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

Basic Bayesian Linear Regression prediction with PyMC3

linear-regression

pymc3