多层次模型+pymc3.glm

Question

我需要使用PyMC3拟合一个多级线性模型，我真的很喜欢glm-api，因为它提供了简洁性。我想问一下是否可以以及如何做到这一点。我发现的这篇博客文章提到：

glm((在分层模型中表现不佳，但

所以我有点怀疑这是否真的可以做到，但它是几年前写的，所以这可能已经改变了。举个例子，下面是我想使用glmAPI重写的模型

import numpy as np
import pymc3 as pm

def generate_data():
n, beta_0, beta_1, sd_eps = 100, 1.2, 0.6, 0.2
b_group = np.array([0.05, 0.14, -0.23])
x = np.random.randn(n)
group_index = np.random.choice([0, 1, 2], n)
y = 1.2 + 0.6 * x + sd_eps * np.random.randn(n) + b_group[group_index]
return x, group_index, y

if __name__ == '__main__':
x, group_index, y = generate_data()
with pm.Model() as multi_level_model:
sd_b_group = pm.HalfNormal("sd_b_group", sigma=100)
b_group = pm.Normal("b_group", mu=0, sigma=sd_b_group, shape=3)
beta_0 = pm.Normal("beta_0", mu=0, sigma=100)
beta_1 = pm.Normal("beta_1", mu=0, sigma=100)
sd_eps = pm.HalfNormal("sd_eps", sigma=100)
pm.Normal("y", beta_0 + beta_1 * x + b_group[group_index],
sigma=sd_eps, observed=y)

我想它看起来像这样：

with pm.Model():
mu_b_group = pm.Normal("mu_b_group", mu=0, sigma=100)
sd_b_group = pm.HalfNormal("sd_b_group", sigma=100)
b_group = pm.Normal("b_group", mu=mu_b_group, sigma=sd_b_group, shape=3)
pm.glm.GLM.from_formula(formula="y ~ 1 + x",
vars={"Intercept": b_group[group_index]},
data={"y": y, "x": x})

然而，当试图堆叠系数时

TypeError: Join() can only join tensors with the same number of dimensions.

Answer 1

经过一些实验，我想出了这个(不理想，但至少有点有用(解决方案：

with pm.Model():
sd_b_group = pm.HalfNormal("sd_b_group", sigma=100)
b_group = pm.Normal("b_group", mu=0, sigma=sd_b_group, shape=3)
lm = pm.glm.LinearComponent.from_formula(formula="y ~ 1 + x",
data={"y": y, "x": x})
sd_eps = pm.HalfNormal("sd_eps", sigma=100)
likelihood = pm.Normal("y", mu=lm.y_est + b_group[group_index], 
sigma=sd_eps, observed=y)        

还在github上创建了一个示例以供将来参考。