在 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 class。）

如有任何帮助，我们将不胜感激。

编辑：添加一些模拟数据

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：

似乎有几种不同的方法可以解决这个问题，尽管 none 完全令人满意：

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

无论如何，开发人员 advice/input 将不胜感激。

Answer 1

少量更改（注意我将 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 )

Answer 2

事实证明解决方案很简单：似乎任何分布（如 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 3

在您的问题中，您已声明 yhat。您可以避免它并将等式传递给 Bernoulli.

的 logit_p 参数

注意 - 您可以传递 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)

在 pymc3 中创建三级逻辑回归模型

Creating a three-level logistic regression model in pymc3

bayesian

pymc