python - PyMC3 将随机协方差矩阵传递给 pm.MvNormal()

Question

我尝试使用 PyMC3 将简单的 2D 高斯模型拟合到观察到的数据。

import numpy as np
import pymc3 as pm

n = 10000;
np.random.seed(0)
X = np.random.multivariate_normal([0,0], [[1,0],[0,1]], n);

with pm.Model() as model:
    # PRIORS
    mu = [pm.Uniform('mux', lower=-1, upper=1), 
          pm.Uniform('muy', lower=-1, upper=1)]
    cov = np.array([[pm.Uniform('a11', lower=0.1, upper=2), 0],
                    [0, pm.Uniform('a22', lower=0.1, upper=2)]])

    # LIKELIHOOD
    likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)

with model:
    trace = pm.sample(draws=1000, chains=2, tune=1000)

虽然我可以通过将 sd 传递给 pm.Normal 来在一维中完成此操作，但在将协方差矩阵传递给 pm.MvNormal< 时遇到了一些麻烦.

我哪里出错了？

Answer 1

PyMC3 分布对象不是简单的数字对象或 numpy 数组。相反，它们是 theano 计算图中的节点，通常需要来自 pymc3.math 的操作。或theano.tensor来操纵他们。此外，将 PyMC3 对象放入 numpy 数组中是不必要的，因为它们已经是多维的。

原始模型

根据代码的意图，工作版本将类似于

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

N = 10000
np.random.seed(0)
X = np.random.multivariate_normal(np.zeros(2), np.eye(2), size=N)

with pm.Model() as model:
    # use `shape` argument to define tensor dimensions
    mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

    # diagonal values on covariance matrix
    a = pm.Uniform('a', lower=0.1, upper=2, shape=2)

    # convert vector to a 2x2 matrix with `a` on the diagonal
    cov = tt.diag(a)

    likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)

替代模型

我认为您提供的示例只是一个传达问题的玩具。但为了以防万一，我会提到，当将协方差矩阵限制为对角矩阵时，使用多元正态(参数之间的协方差建模)的主要优点就会丢失。此外，协方差矩阵的先验理论已经很成熟，因此值得花时间考虑现有的解决方案。特别是有a PyMC3 example using the LKJ prior for covariance matrices .

以下是该示例在此上下文中的简单应用:

with pm.Model() as model_lkj:
    # use `shape` argument to define tensor dimensions
    mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

    # LKJ prior for covariance matrix (see example)
    packed_L = pm.LKJCholeskyCov('packed_L', n=2,
                                 eta=2., sd_dist=pm.HalfCauchy.dist(2.5))
    # convert to (2,2)
    L = pm.expand_packed_triangular(2, packed_L)

    likelihood = pm.MvNormal('likelihood', mu=mu, chol=L, observed=X)