使用 PYMC3 进行分层线性回归的多层次

Question

我正在尝试使用 PYMC3 建立分层线性回归模型。在我的特殊情况下，我想看看邮政编码是否为其他功能提供了有意义的结构。假设我使用以下模拟数据：

import pandas as pd
import numpy as np
import pymc3 as pm

data = pd.DataFrame({"postalcode": np.floor(np.random.uniform(low=10, high=99, size=1000)),
                 "x": np.random.normal(size=1000),
                 "y": np.random.normal(size=1000)})
data["postalcode"] = data["postalcode"].astype(int)

我生成从 10 到 99 的邮政编码，以及一个正态分布的特征 x 和一个目标值 y。现在我为邮政编码级别 1 和级别 2 设置索引：

def create_pc_index(level):
    pc = data["postalcode"].astype(str).str[0:level]
    unique_pc = pc.unique()
    pc_dict = dict(zip(unique_pc, range(0, len(unique_pc))))
    return pc_dict, pc.apply(lambda x: pc_dict[x]).values

pc1_dict, pc1_index = create_pc_index(1)
pc2_dict, pc2_index = create_pc_index(2) 

使用邮政编码的第一位作为分层属性可以正常工作：

number_of_samples = 1000

x = data["x"]
y = data["y"]

with pm.Model() as model:
    sigma = pm.HalfCauchy('sigma', beta=10, testval=0.5, shape=1)
    mu_i = pm.Normal("mu_i", 5, sd=25, shape=1)
    intercept = pm.Normal('Intercept', mu_i, sd=1, shape=len(pc1_dict))

    mu_s = pm.Normal("mu_x", 0, sd=3, shape=1)
    x_coeffs = pm.Normal("x", mu_s, 1, shape=len(pc1_dict))

    mean = intercept[pc1_index] + x_coeffs[pc1_index] * x

    likelihood_mean = pm.Deterministic("mean", mean)
    likelihood = pm.Normal('y', mu=likelihood_mean, sd=sigma, observed=y)

    trace = pm.sample(number_of_samples)
    burned_trace = trace[number_of_samples/2:]

但是，如果我想在我的层次结构中添加第二个级别（在这种情况下仅在截距上，暂时忽略 x），我会遇到形状问题

with pm.Model() as model:
    sigma = pm.HalfCauchy('sigma', beta=10, testval=0.5, shape=1)
    mu_i_level_1 = pm.Normal("mu_i", 0, sd=25, shape=1)
    mu_i_level_2 = pm.Normal("mu_i_level_2", mu_i_level_1, sd=1, shape=len(pc1_dict))
    intercept = pm.Normal('Intercept', mu_i_level_2[pc1_index], sd=1, shape=len(pc2_dict))

    mu_s = pm.Normal("mu_x", 0, sd=3, shape=1)
    x_coeffs = pm.Normal("x", mu_s, 1, shape=len(pc1_dict))

    mean = intercept[pc2_index] + x_coeffs[pc1_index] * x

    likelihood_mean = pm.Deterministic("mean", mean)
    likelihood = pm.Normal('y', mu=likelihood_mean, sd=sigma, observed=y)

    trace = pm.sample(number_of_samples)
    burned_trace = trace[number_of_samples/2:]

错误信息是：

operands could not be broadcast together with shapes (89,) (1000,) 

如何正确建模回归中的多个级别？这只是形状尺寸正确的问题，还是我有更根本的错误？

提前致谢！

Answer 1

我认为intercept 的形状不能是len(pc2_dict)，而是len(pc1_dict) 的mu。矛盾就在这里：

intercept = pm.Normal('Intercept', mu_i_level_2[pc1_index], sd=1, shape=len(pc2_dict))