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

问题描述

我已经尝试通过使用PyMC3将简单的2D高斯模型拟合到观测数据.

 将numpy导入为np导入pymc3作为pmn = 10000；np.random.seed(0)X = np.random.multivariate_normal([0,0]，[[1,0]，[0,1]]，n);使用pm.Model()作为模型:#优先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)]])#喜欢可能性= pm.MvNormal('可能性'，mu = mu，cov = cov，观察到的= X)与模型:跟踪= pm.sample(绘制= 1000，链= 2，调= 1000) 

虽然我可以通过将 sd 传递给 pm.Normal 在1D中做到这一点，但在将协方差矩阵传递给 pm.MvNormal .

我要去哪里错了?

解决方案

PyMC3分发对象不是简单的数字对象或numpy数组.相反，它们是theano计算图中的节点，通常需要来自 pymc3.math的操作 或 theano.tensor 进行操作.此外，由于PyMC3对象已经是多维的，因此不必将它们放置在numpy数组中.

原始模型

按照您的代码意图，工作版本会像这样

 将numpy导入为np导入pymc3作为pm将theano.tensor导入为ttN = 10000np.random.seed(0)X = np.random.multivariate_normal(np.zeros(2)，np.eye(2)，size = N)使用pm.Model()作为模型:#使用`shape`参数定义张量尺寸mu = pm.Uniform('mu'，lower = -1，upper = 1，shape = 2)#协方差矩阵上的对角线值a = pm.Uniform('a'，lower = 0.1，upper = 2，shape = 2)#将向量转换为对角线为a的2x2矩阵cov = tt.diag(a)可能性= pm.MvNormal('可能性'，mu = mu，cov = cov，观察到的= X) 

替代模型

我假设您提供的示例只是一个用来传达问题的玩具.但是以防万一，我要提到的是，当将协方差矩阵限制为对角线时，将失去使用多元法线(参数之间的建模协方差)的主要优势.此外，协方差矩阵的先验理论已经得到了很好的发展，因此考虑现有的解决方案值得花时间.特别是，有使用LKJ优先用于协方差矩阵的PyMC3示例.>

在这种情况下，以下是该示例的简单应用程序:

 ，其中pm.Model()为model_lkj:#使用`shape`参数定义张量尺寸mu = pm.Uniform('mu'，lower = -1，upper = 1，shape = 2)#LKJ先于协方差矩阵(请参见示例)packed_L = pm.LKJCholeskyCov('packed_L'，n = 2，eta = 2.，sd_dist = pm.HalfCauchy.dist(2.5))#转换为(2,2)L = pm.expand_packed_triangular(2，packed_L)可能性= pm.MvNormal('可能性'，mu = mu，chol = L，观察= X) 

I've tried to fit a simple 2D gaussian model to observed data by using PyMC3.

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)

while I can do this in 1D by passing the sd to pm.Normal, I have some trouble in passing the covariance matrix to pm.MvNormal.

Where am I going wrong?

解决方案

PyMC3 distribution objects are not simple numeric objects or numpy arrays. Instead, they are nodes in a theano computation graph and often require operations from either pymc3.math or theano.tensor to manipulate them. Moreover, placing PyMC3 objects in numpy arrays is unnecessary since they already are multidimensional.

Original Model

Keeping with the intent of your code, a working version would go something like

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)

Alternative Model

I assume the example you provided is just a toy to communicate the problem. But just in case, I'll mention that the main advantage of using a multivariate normal (modeling covariance between parameters) is lost when restricting the covariance matrix to be diagonal. Furthermore, the theory of priors for covariance matrices is well-developed, so it's worth one's time considering existing solutions. In particular, there is a PyMC3 example using the LKJ prior for covariance matrices.

Here's a simple application of that example in this context:

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)

这篇关于PyMC3将随机协方差矩阵传递给pm.MvNormal()的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！