Создание трехуровневой модели логистической регрессии в pymc3

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

Здесь мой код:

with pm.Model() as model:
    # Hyperpriors
    top_level_tau = pm.HalfNormal('top_level_tau', sd=100.)
    mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)    

    # Priors
    top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top)
    mid_level = [pm.Normal('mid_level_{}'.format(j),
                           mu=top_level[mid_to_top_idx[j]],
                           tau=mid_level_tau)
                 for j in range(k_mid)]

    intercept = pm.Normal('intercept', mu=0., sd=100.)

    # Model prediction
    yhat = pm.invlogit(mid_level[mid_to_bot_idx] + intercept)

    # Likelihood
    yact = pm.Bernoulli('yact', p=yhat, observed=y)

Я получаю сообщение об ошибке "only integer arrays with one element can be converted to an index" (в строке 16), которое, как мне кажется, связано с тем, что переменная mid_level является списком, а не соответствующим контейнером pymc. (Я не вижу класс Container в исходном коде pymc3.)

Любая помощь будет оценена.

Изменить: добавление некорректных данных

y = np.array([0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0])
mid_to_bot_idx = np.array([0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 3, 2, 3, 3, 3, 3, 2, 2, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 2, 2, 3, 2, 2, 3, 3, 3, 3, 2, 2, 2, 3, 2, 3, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 2])
mid_to_top_idx = np.array([0, 0, 1, 1])
k_top = 2
k_mid = 4

Изменить # 2:

Кажется, существует несколько способов решения этой проблемы, хотя ни один из них не является полностью удовлетворительным:

1) Можно пересмотреть модель следующим образом:

with pm.Model() as model:
    # Hyperpriors
    top_level_tau = pm.HalfNormal('top_level_tau', sd=100.)
    mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)    

    # Priors
    top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top)
    mid_level = pm.Normal('mid_level', mu=0., tau=mid_level_tau, shape=k_top)
    intercept = pm.Normal('intercept', mu=0., sd=100.)

    # Model prediction
    yhat = pm.invlogit(top_level[top_to_bot_idx] + mid_level[mid_to_bot_idx] + intercept)

    # Likelihood
    yact = pm.Bernoulli('yact', p=yhat, observed=y)

Кажется, что это работает, хотя я не могу понять, как его распространить на случай, когда дисперсия среднего уровня не является постоянной для всех групп среднего уровня.

2) Можно скомпоновать коэффициенты среднего уровня в тензор Anano с использованием theano.tensor.stack: i.e.,

import theano.tensor as tt
mid_level = tt.stack([pm.Normal('mid_level_{}'.format(j),
                           mu=top_level[mid_to_top_idx[j]],
                           tau=mid_level_tau)
                 for j in range(k_mid)])

Но это кажется очень медленным в моем фактическом наборе данных (30 тыс. наблюдений), и это делает невозможным построение графика (каждый из коэффициентов среднего уровня получает свою собственную трассировку, используя pm.traceplot).

В любом случае, некоторые советы/вклад от разработчиков будут оценены.

Ответы

Ответ 1

Оказывается, решение было простым: кажется, что любое распределение (например, pm.Normal) может принимать вектор значений в качестве аргумента, поэтому замена этой строки

mid_level = [pm.Normal('mid_level_{}'.format(j),
                       mu=top_level[mid_to_top_idx[j]],
                       tau=mid_level_tau)
             for j in range(k_mid)]

с этим

mid_level = pm.Normal('mid_level',
                       mu=top_level[mid_to_top_idx],
                       tau=mid_level_tau,
                       shape=k_mid)

работает. Тот же метод также может использоваться для указания отдельных стандартных отклонений для каждой из групп среднего уровня.

EDIT: исправлена опечатка

Ответ 2

Немногие изменения (обратите внимание, что я изменил yhat на theta):

theta = pm.Deterministic( 'theta', pm.invlogit(sum(mid_level[i] for i in mid_to_bot_idx)+intercept) )
yact = pm.Bernoulli( 'yact', p=theta, observed=y )

Ответ 3

В своем вопросе вы заявили о yhat. Вы можете избежать этого и передать уравнение параметру logit_p Bernoulli.

Примечание. Вы можете передать либо p либо logit_p.

В моем случае использование logit_p ускорит процесс выборки.

Code-

# Likelihood
yact = pm.Bernoulli('yact', logit_p=top_level[top_to_bot_idx] + mid_level[mid_to_bot_idx] + intercept, observed=y)