python - 我们如何在 PyMC3 的分层模型中预测新的看不见的组？

Question

如果我们有一个分层模型，其中来自不同站点的数据作为模型中的不同组，那么我们如何预测新组(我们以前没有见过的新站点)？例如使用以下逻辑回归模型:

from pymc3 import Model, sample, Normal, HalfCauchy,Bernoulli
import theano.tensor as tt
with Model() as varying_slope:

    mu_beta = Normal('mu_beta', mu=0., sd=1e5)
    sigma_beta = HalfCauchy('sigma_beta', 5)
    a = Normal('a', mu=0., sd=1e5)
    betas = Normal('b',mu=mu_beta,sd=sigma_beta,shape=(n_features,n_site))
    y_hat = a + tt.dot(X_shared,betas[:,site_shared])
    y_like = Bernoulli('y_like', logit_p=y_hat, observed=train_y)

拟合此模型后，我们可以使用以下方法对来自特定站点的新数据(即后验预测的样本)进行预测:

site_to_predict = 1
samples = 500

x = tt.matrix('X',dtype='float64')
new_site = tt.vector('new_site',dtype='int32')
n_samples = tt.iscalar('n_samples')
x.tag.test_value = np.empty(shape=(1,X.shape[1]))
new_site.tag.test_value = np.empty(shape=(1,1))

_sample_proba =  approx.sample_node(varying_slope.y_like.distribution.p,
                               size=n_samples,
                               more_replacements={X_shared: x,site_shared:new_site})

sample_proba = theano.function([x,new_site,n_samples], _sample_proba)

pred_test = sample_proba(test_X.reshape(1,-1),np.array(site_to_predict).reshape(-1),samples)

但是如果我们有一个新的未见过的站点，从后验预测分布中采样的正确方法是什么？

Answer 1

我只是从 pymc discourse thread 复制我的答案如果有人偶然遇到这个问题或另一个类似的问题。

首先，请注意您正在使用的居中分层参数化1，它可能会导致拟合时出现分歧和困难。

话虽这么说，您的模型看起来或多或少像 GLM，具有跨功能和站点共享先验随机变量 mu_beta 和 sigma_beta。一旦你获得了这两者的后验分布，你的预测应该类似于

y_hat = a + dot(X_shared, Normal(mu=mu_beta, sigma=sigma_beta))
y_like = Bernoulli('y_like', logit_p=y_hat)

因此，我们的目标是实现这一目标。

我们始终推荐样本外后验预测检查的方式是使用 theano.shared。我将使用一种不同的方法，其灵感来自于作为 pymc4 核心设计理念的函数式 API。我不会讨论 pymc3 和 pymc4 的骨架之间的许多差异，但我开始更多地使用工厂函数来获取模型实例。我没有尝试使用 theano.shared 定义模型内部的内容，而是使用新数据创建一个新模型并从中提取后验预测样本。我最近刚刚在这里发布了相关内容。

这个想法是使用训练数据创建模型并从中采样以获得跟踪。然后，您必须从跟踪中提取与看不见的站点共享的分层部分:mu_beta、sigma_beta 和 a。最后，您使用测试站点的新数据创建一个新模型，并使用保存 mu_beta、sigma_beta 和部分训练轨迹的字典列表从后验预测中进行采样。这是一个独立的示例

import numpy as np
import pymc3 as pm
from theano import tensor as tt
from matplotlib import pyplot as plt


def model_factory(X, y, site_shared, n_site, n_features=None):
    if n_features is None:
        n_features = X.shape[-1]
    with pm.Model() as model:
        mu_beta = pm.Normal('mu_beta', mu=0., sd=1)
        sigma_beta = pm.HalfCauchy('sigma_beta', 5)
        a = pm.Normal('a', mu=0., sd=1)
        b = pm.Normal('b', mu=0, sd=1, shape=(n_features, n_site))
        betas = mu_beta + sigma_beta * b
        y_hat = a + tt.dot(X, betas[:, site_shared])
        pm.Bernoulli('y_like', logit_p=y_hat, observed=y)
    return model


# First I generate some training X data
n_features = 10
ntrain_site = 5
ntrain_obs = 100
ntest_site = 1
ntest_obs = 1
train_X = np.random.randn(ntrain_obs, n_features)
train_site_shared = np.random.randint(ntrain_site, size=ntrain_obs)
new_site_X = np.random.randn(ntest_obs, n_features)
test_site_shared = np.zeros(ntest_obs, dtype=np.int32)
# Now I generate the training and test y data with a sample from the prior
with model_factory(X=train_X,
                   y=np.empty(ntrain_obs, dtype=np.int32),
                   site_shared=train_site_shared,
                   n_site=ntrain_site) as train_y_generator:
    train_Y = pm.sample_prior_predictive(1, vars=['y_like'])['y_like'][0]
with model_factory(X=new_site_X,
                   y=np.empty(ntest_obs, dtype=np.int32),
                   site_shared=test_site_shared,
                   n_site=ntest_site) as test_y_generator:
    new_site_Y = pm.sample_prior_predictive(1, vars=['y_like'])['y_like'][0]

# The previous part is just to get some toy data to fit
# Now comes the important parts. First training
with model_factory(X=train_X,
                   y=train_Y,
                   site_shared=train_site_shared,
                   n_site=ntrain_site) as train_model:
    train_trace = pm.sample()

# Second comes the hold out data posterior predictive
with model_factory(X=new_site_X,
                   y=new_site_Y,
                   site_shared=test_site_shared,
                   n_site=ntrain_site) as test_model:
    # We first have to extract the learnt global effect from the train_trace
    df = pm.trace_to_dataframe(train_trace,
                               varnames=['mu_beta', 'sigma_beta', 'a'],
                               include_transformed=True)
    # We have to supply the samples kwarg because it cannot be inferred if the
    # input trace is not a MultiTrace instance
    ppc = pm.sample_posterior_predictive(trace=df.to_dict('records'),
                                         samples=len(df))
    plt.figure()
    plt.hist(ppc['y_like'], 30)
    plt.axvline(new_site_Y, linestyle='--', color='r')

我得到的后验预测如下所示:

当然，我不知 Prop 体应该把什么样的数据作为你的X_shared、site_shared或train_y，所以我只是在代码开头编造了一些无意义的玩具数据，你应该用你的实际数据替换它.