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