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

Question

Я пытаюсь создать трехуровневую модель логистической регрессии в 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.)

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

Edit: Добавление некоторых фиктивные данные

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 

Edit # 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) Можно обернуть коэффициенты среднего уровня в тензор Theano использование theano.tensor.stack: т.е.

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)]) 

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

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

Answer 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: Исправлена ​​опечатка

Answer 2

Некоторые изменения (обратите внимание, я изменил yhat на тета):

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)