在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容器。 （我沒有在pymc3源代碼中看到Container類。）

任何幫助，將不勝感激。

編輯：添加一些模擬數據

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），可以使用theano.tensor.stack包裹中級係數爲Theano張量：即

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) 

作品。也可以使用相同的方法來爲每個中等水平的羣體指定單獨的標準偏差。

編輯：固定錯字

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