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

Question

У меня возникает общая проблема. Мне интересно, с кем можно помочь. Я часто хотел бы использовать 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.)

Answer 1

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

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

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 работает очень медленно. Это занимает столько же времени, сколько и фактический вывод. Я подозреваю, что это связано с некоторой неэффективностью, связанной с использованием анано.

EDIT: ссылка, первоначально включенная, кажется, мертва.