Создание прогнозов на основе предполагаемых параметров в 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]), которая просто пропускала бы новые данные через граф вычислений theano.

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

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


person santon    schedule 21.10.2015    source источник
comment
Вы говорите о выборке из апостериорного прогнозирующего распределения, что, по-видимому, вы делаете правильно. не совсем уверен, что вы имеете в виду, основываясь на новых данных. Вы говорите об использовании апостериорных значений из этого анализа в качестве априорных для вывода на основе дополнительных данных?   -  person Chris Fonnesbeck    schedule 21.10.2015
comment
@ChrisFonnesbeck Это то, что мне также было бы интересно, поскольку апостериорные данные, которые мы получаем, находятся в форме трассировки, и мы не можем использовать их для указания априорных значений в синтаксисе примера.   -  person recluze    schedule 21.10.2015
comment
twiecki на странице gitter pymc3 указал мне на эту страницу, которая, похоже, касается проблема, с которой я столкнулся. Придется потратить некоторое время, чтобы понять, что было сделано, но это выглядит многообещающе.   -  person santon    schedule 21.10.2015
comment
При достаточном количестве данных апериорные зубы обычно (многомерны) нормальны. Простым подходом было бы извлечь среднее значение и стандартное отклонение апостериорного значения и использовать это для параметризации нормальных априорных значений для последующего анализа.   -  person Chris Fonnesbeck    schedule 22.10.2015


Ответы (2)


Примечание. Эта функция теперь включена в основной код как метод 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.

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

person santon    schedule 22.10.2015

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

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

person Ashok Rayal    schedule 26.04.2019