Формирование прогнозов из предполагаемых параметров в pymc3

Я сталкиваюсь с общей проблемой, о которой я спрашиваю, может ли кто-то помочь. Я часто хотел бы использовать pymc3 в двух режимах: обучение (то есть фактическое выполнение вывода по параметрам) и оценка (т.е. Использование предполагаемых параметров для генерации предсказаний).

В общем, мне хотелось бы, чтобы апостериор по предсказаниям, а не только по подсчетам (это часть преимуществ байесовской рамки, нет?). Когда ваши данные обучения фиксированы, это обычно выполняется путем добавления моделируемой переменной подобной формы к наблюдаемой переменной. Например,

from pymc3 import *

with basic_model:

    # Priors for unknown model parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
    Y_sim = Normal('Y_sim', mu=mu, sd=sigma, shape=len(X1))

    start = find_MAP()
    step = NUTS(scaling=start)
    trace = sample(2000, step, start=start)

Но что, если мои данные меняются? Скажем, я хочу генерировать прогнозы на основе новых данных, но без повторения вывода снова и снова. В идеале у меня была бы функция типа predict_posterior(X1_new, X2_new, 'Y_sim', trace=trace) или даже predict_point(X1_new, X2_new, 'Y_sim', vals=trace[-1]), которая просто запускала бы новые данные с помощью графика вычислений anano.

Я полагаю, что часть моего вопроса связана с тем, как pymc3 реализует диаграмму вычисления anano. Я заметил, что функция model.Y_sim.eval кажется похожей на то, что я хочу, но она требует Y_sim в качестве ввода и, кажется, просто возвращает все, что вы ей даете.

Я предполагаю, что этот процесс чрезвычайно распространен, но я не могу найти способ сделать это. Любая помощь приветствуется. (Обратите внимание, что у меня есть взломать это в pymc2, это сложнее в pymc3 из-за theano.)

Ответ 1

Примечание: эта функциональность теперь включена в основной код как метод pymc.sample_ppc. Проверьте документы для получения дополнительной информации.

Основываясь на этой ссылке (мертвой по состоянию на июль 2017 года), присланной мне twiecki, есть несколько хитростей, чтобы решить мою проблему. Первый - поместить данные обучения в общую переменную theano. Это позволяет нам позже изменить данные, не облажаясь с графиком анано вычислений.

X1_shared = theano.shared(X1)
X2_shared = theano.shared(X2)

Затем создайте модель и выполните вывод как обычно, но с использованием общих переменных.

with basic_model:

    # Priors for unknown model parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1_shared + beta[1]*X2_shared

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

    start = find_MAP()
    step = NUTS(scaling=start)
    trace = sample(2000, step, start=start)

Наконец, есть функция в стадии разработки (которая, вероятно, в конечном итоге будет добавлена в pymc3), которая позволит прогнозировать исходные данные для новых данных.

from collections import defaultdict

def run_ppc(trace, samples=100, model=None):
    """Generate Posterior Predictive samples from a model given a trace.
    """
    if model is None:
         model = pm.modelcontext(model)

    ppc = defaultdict(list)
    for idx in np.random.randint(0, len(trace), samples):
        param = trace[idx]
        for obs in model.observed_RVs:
            ppc[obs.name].append(obs.distribution.random(point=param))

    return ppc

Затем передайте новые данные, для которых вы хотите выполнить прогнозы:

X1_shared.set_value(X1_new)
X2_shared.set_value(X2_new)

Наконец, вы можете сгенерировать апостериорные прогностические выборки для новых данных.

ppc = run_ppc(trace, model=model, samples=200)

Переменная ppc - это словарь с ключами для каждой наблюдаемой переменной в модели. Таким образом, в этом случае ppc['Y_obs'] будет содержать список массивов, каждый из которых генерируется с использованием одного набора параметров из трассировки.

Обратите внимание, что вы даже можете изменить параметры, извлеченные из трассировки. Например, у меня была модель, использующая переменную GaussianRandomWalk и я хотел генерировать прогнозы на будущее. В то время как вы можете позволить pymc3 выполнять выборку в будущем (т.е. Позволить переменной случайного блуждания расходиться), я просто хотел использовать фиксированное значение коэффициента, соответствующего последнему выведенному значению. Эта логика может быть реализована в функции run_ppc.

Стоит также отметить, что функция run_ppc работает крайне медленно. Это занимает примерно столько же времени, сколько и фактический вывод. Я подозреваю, что это связано с некоторой неэффективностью, связанной с использованием theano.

РЕДАКТИРОВАТЬ: Ссылка изначально включена, кажется, мертв.

Ответ 2

Выше ответ от @santon является правильным. Я просто добавляю к этому.

Теперь вам не нужно писать свой собственный метод run_ppc. pymc3 предоставляет метод sample_posterior_predictive который делает то же самое.