1

如果我们有一个分层模型，将来自不同站点的数据作为模型中的不同组，我们如何预测新组（我们以前从未见过的新站点）？例如使用以下逻辑回归模型：

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)

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

4

1 回答 1

3

如果有人偶然遇到这个问题或其他类似的问题，我只是从pymc 讨论线程中复制我的答案。

首先，请注意您正在使用的居中分层参数化 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's。我将使用不同的方法，灵感来自作为 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')

我得到的后验预测如下所示： 在此处输入图像描述

当然，我不知道你的X_shared、site_shared还是train_y具体放什么数据，所以我只是在代码的开头编了一些废话玩具数据，你应该用你的实际数据替换它。

于 2019-01-22T22:03:20.383 回答