具有 MVN 先验的分层变化效应模型

Question

我想做什么

我已经在 pymc 中处理过多元先验（我使用的是 4.0.1），但我无法在分层模型中使用它们。在我的示例中，我使用两个协变量

x1

、

x2

和结果

y

来建模回归问题。数据中有两个定义层次结构级别的分类特征：

dim1

表示较高级别类别，

dim2

表示较低级别类别。

该模型由截距和斜率组成，分别对应

x1

和

x2

。有一个多元正态超先验

B_1

，它按

dim1

的顺序定义变化的截距和斜率 [intercept,lope_x1,slope_x2] 。然后，这些信息通知较低级别上的先前

B_2

，其中术语也因

dim2

而变化。

问题

尝试从模型中采样时，出现以下错误：“ValueError：值的维度无效：3”。
我将其解释为

B_1

和

B_2

的暗淡中的“3”是错误的，但是当手动评估模型变量时，一切似乎都很好：

B_2.eval().shape

正确显示 rv 的形状为 (6,2,3)。

mu

中的线性模型也可以毫无问题地进行评估。

检查错误回溯时提到了这个问题：

/local_disk0/.ephemeral_nfs/cluster_libraries/python/lib/python3.9/site-packages/pymc/distributions/multivariate.py in quaddist_parse(value, mu, cov, mat_type)
    133     """Compute (x - mu).T @ Sigma^-1 @ (x - mu) and the logdet of Sigma."""
    134     if value.ndim > 2 or value.ndim == 0:
--> 135         raise ValueError("Invalid dimension for value: %s" % value.ndim)
    136     if value.ndim == 1:
    137         onedim = True

这意味着维度不能大于 2。这看起来很奇怪，因为维度超过 2 的多元分布非常常见。这让我相信我指定 MVN 先验的方式是错误的，但经过多次尝试后我陷入了困境。 在这里指定

B_1

和

B_2

的正确方法是什么？

型号代码

import numpy as np
import pymc as pm
with pm.Model(coords={'dim1':np.arange(2), 'dim2':np.arange(6)}) as mdl:
    d1 = pm.MutableData('d1', X['d1'].values, dims='obs_id')
    d2 = pm.MutableData('d2', X['d2'].values, dims='obs_id')
    # upper level
    sd_1 = pm.HalfNormal.dist(1,shape=3)
    chol_1, _, _ = pm.LKJCholeskyCov('chol_1', n=3, eta=1,
        sd_dist=sd_1, compute_corr=True)
    B_1 = pm.MvNormal('B_1', mu=[5,0,0], chol=chol_1, 
                        dims=('dim1','3')
                        )
    # lower level
    sd_2 = pm.HalfNormal.dist(1,shape=3)
    chol_2, _, _ = pm.LKJCholeskyCov('chol_2', n=3, eta=1,
        sd_dist=sd_2, compute_corr=True)
    B_2 = pm.MvNormal('B_2', mu=B_1, chol=chol_2, 
                        dims=('dim2','dim1','3')
                        )
    # regular robust regression, mu defined as linear model of covariates
    sigma = pm.HalfNormal('sigma',2)
    mu = pm.Deterministic(
        'mu', 
        B_2[d2,d1,0] + # intercept
        B_2[d2,d1,1] * X['x1'].values + # slope x1
        B_2[d2,d1,2] * X['x2'].values, # slope x2
        dims='obs_id'
    )
    outcome = pm.StudentT('outcome', nu=3, mu=mu, sigma=sigma, observed=y, dims="obs_id")
    trace = pm.sample(draws=2000, tune=1000, target_accept=0.95, random_seed=0)

数据生成代码

import numpy as np
import pandas as pd
from scipy import stats
import pymc as pm

rng = np.random.default_rng(0)
N = 1000
# generate covariates and categories
X = pd.DataFrame(
    {
        'x1': stats.halfnorm(loc=0,scale=3).rvs(N),
        'x2': stats.norm(loc=0,scale=2).rvs(N),
        'd1': rng.choice([0,1],size=N, p=[0.6,0.4]),
        'd2': rng.choice(np.arange(6),size=N, p=[0.1,0.2,0.1,0.3,0.1,0.2]),
    }
)
# means of the parameter distributions
intercept = np.array([
    [5,5,6,7,6,5],
    [4,4,5,6,6,4]
])
slope1 = np.array([
    [0,0.7,0.3,-0.2,-1,0],
    [0.5,1,-1,0.6,-0.2,0.3]
])
slope2 = np.array([
    [0,0.7,0.3,-0.2,-1,0],
    [0.5,1,-1,0.6,-0.2,0.3]
])*1.5
# generate some random covariance matrices
corrs = []
for _ in np.arange(6):
    _,corr,_ = pm.LKJCholeskyCov.dist(eta=1,n=3,sd_dist=pm.HalfNormal.dist(1,shape=3), compute_corr=True)
    corrs.append(corr.eval())

# generate outcome 
y = np.zeros(N)
for d1 in [0,1]:
    for d2 in np.arange(6):
        ind = (X['d1']==d1)&(X['d2']==d2)
        mv = stats.multivariate_normal(mean=[intercept[d1,d2], slope1[d1,d2],slope2[d1,d2]],cov=corrs[d2]).rvs(1)
        y[ind] = mv[0] + X.loc[ind,'x1']*mv[1] + X.loc[ind,'x2']*mv[2] + rng.normal(loc=0,scale=1,size=ind.sum())

Answer 1

好吧，在用尽所有明显的选项后，我发现 PyMC 目前不支持此功能，但如果支持的话，该方法本身就很好。希望我可以通过这种方式节省别人的时间。