多层模型+ pymc3.glm

问题描述

我需要使用PyMC3拟合多级线性模型，并且我非常喜欢glm api，因为它提供的简洁性。我想问一下是否以及如何做到这一点。我发现的blog post提到：

glm（）不适用于分层模型

所以我有点怀疑这实际上可以完成，但是它是几年前写的，所以可能已经改变了。举个例子，下面是我想使用glm api重写的模型

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,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",sigma=100)
        beta_1 = pm.normal("beta_1",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",sigma=100)
    sd_b_group = pm.Halfnormal("sd_b_group",sigma=100)

    b_group = pm.normal("b_group",mu=mu_b_group,shape=3)

    pm.glm.GLM.from_formula(formula="y ~ 1 + x",vars={"Intercept": b_group[group_index]},data={"y": y,"x": x})

但是在尝试堆叠系数时会产生错误here

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

解决方法

经过一些试验，我想出了这个（不理想，但至少有用的）解决方案：

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)