具有分离均值、维度不匹配错误的多元正态模型

Model multivariate normal with separate means, dimension mismatch error

我想用单独的均值和协方差矩阵对多变量(示例中的双变量)法线建模,但不知何故我无法获得匹配的维度。模型比较简单,下面是数据生成过程:

import numpy as np
import pymc3 as pm
import theano.tensor as tt

np.random.seed(123)
N = 200
# x1 and x2 are independent variables
x1 = np.random.randn(N)
x2 = np.random.randn(N)
α1, β1, α2, β2 = np.random.randn(4)

# the means of the joint distribution
μ1 = α1 + β1 * x1
μ2 = α2 + β2 * x2

# the covariance of the joint distribution
σ = 1, 2
Ω = np.array([[1, 0.5], [0.5, 1]])
Σ = np.diag(σ) @ Ω @ np.diag(σ)

# the error term 
ϵ = np.random.multivariate_normal(mean=(0, 0), cov=Σ, size=N)

# the dependent variable as bivariate normal
y = np.stack((μ1, μ2)).T + ϵ

这是我的模型:

with pm.Model() as joint_model:
    a1 = pm.Normal('a1', 0, 1)
    b1 = pm.Normal('b1', 0, 1)
    a2 = pm.Normal('a2', 0, 1)
    b2 = pm.Normal('b2', 0, 1)

    mu1 = a1 + b1 * x1
    mu2 = a2 + b2 * x2

    sigma = pm.HalfCauchy('sigma', 1, shape=2)

    # build the covariance matrix
    C_triu = pm.LKJCorr('omega', n=2, p=2)
    C = tt.fill_diagonal(C_triu[np.zeros((2, 2), dtype=np.int64)], 1)
    sigma_diag = tt.nlinalg.diag(sigma)
    cov = tt.nlinalg.matrix_dot(sigma_diag, C, sigma_diag)

    joint_obs = pm.MvNormal('joint', mu=(mu1, mu2), cov=cov, observed=y)

错误信息是ValueError: Input dimension mis-match. (input[0].shape[0] = 200, input[1].shape[0] = 2)。 很明显,我在指定多元正态分布时犯了一个错误,mucov 和观察到的数据之间存在(可能不止一个)维度不匹配,我只是想不通在哪里以及如何修复它。这是(很长的)错误消息,上面的代码应该能够复制错误:

Traceback (most recent call last):

  File "<ipython-input-4-b3658e712eaf>", line 18, in <module>
    joint_obs = pm.MvNormal('joint', mu=(mu1, mu2), cov=cov, observed=y)

  File "/Users/oma/anaconda/lib/python3.6/site-packages/pymc3/distributions/distribution.py", line 39, in __new__
    return model.Var(name, dist, data, total_size)

  File "/Users/oma/anaconda/lib/python3.6/site-packages/pymc3/model.py", line 545, in Var
    total_size=total_size, model=self)

  File "/Users/oma/anaconda/lib/python3.6/site-packages/pymc3/model.py", line 970, in __init__
    self.logp_elemwiset = distribution.logp(data)

  File "/Users/oma/anaconda/lib/python3.6/site-packages/pymc3/distributions/multivariate.py", line 200, in logp
    delta = value - mu

  File "/Users/oma/anaconda/lib/python3.6/site-packages/theano/tensor/var.py", line 147, in __sub__
    return theano.tensor.basic.sub(self, other)

  File "/Users/oma/anaconda/lib/python3.6/site-packages/theano/gof/op.py", line 674, in __call__
    required = thunk()

  File "/Users/oma/anaconda/lib/python3.6/site-packages/theano/gof/op.py", line 843, in rval
    fill_storage()

  File "/Users/oma/anaconda/lib/python3.6/site-packages/theano/gof/cc.py", line 1698, in __call__
    reraise(exc_type, exc_value, exc_trace)

  File "/Users/oma/anaconda/lib/python3.6/site-packages/six.py", line 686, in reraise
    raise value

ValueError: Input dimension mis-match. (input[0].shape[0] = 200, input[1].shape[0] = 2)

问题已在 discourse.pymc.io 上得到解答。

pm.MvNormalmu 的规格应为 mu = tt.stack([mu1, mu2]).T