Я просто копирую свой ответ из нити дискурса 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.Я буду использовать другой подход, основанный на функциональном API, который является основной идеей дизайна для pymc4.Существует много различий между pymc3 и каркасом pymc4, но я начал использовать больше, это фабричные функции для получения экземпляров Model.Вместо того, чтобы пытаться определять вещи внутри модели с помощью 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, поэтому я просто собрал некоторые бессмысленные игрушечные данные в начале кода, вы должны заменить их фактическими данными.