在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观测值），并且它使得绘图不方便（每个mid_level系数使用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

几点变化（注意我将其更改为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 )

Answer 3

在你的问题中，你已经宣布了yhat。你可以避免它并将公式传递给logit_p的Bernoulli参数。

注意 - 您可以传递p或logit_p。

在我的情况下使用logit_p加快我的采样过程。

码-

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