Construct matrix with random variables as elements

I am trying to create an observed variable in PyMC3 which has a multivariate normal distribution. I want the covariance matrix whose elements are other random variables. As an example, consider the following code:

import pymc3 as pm

with pm.Model() as model:
  a = pm.Normal('a', mu=0, sigma=10)
  b = pm.Normal('a', mu=0, sigma=10)
  c = pm.Normal('a', mu=0, sigma=10)

  # I recognize that the matrix is not positive definite
  # with this parameterization. This is just a toy example.
  # The main point is that I want the elements of the matrix
  # to be random variables.
  cov = [[a, b]
         [b, c]]

  data = pm.MvNormal('data', mu=[0, 0],
                       cov=cov, observed=obs)

This does not work. Neither does using np.array(cov). I imagine the solution is to use Theano tensors somehow. I am unable to figure out how to use them.

I’d appreciate any help with this. Thanks.

You can use

import theano.tensor as tt
tt.concatenate([a, b, b, c]).reshape((2, 2))

Or alternatively tt.stack.
If you want priors on covariance matrices, you probably want to look into pm.LKJCholeskyCov instead.

1 Like

Using tt.stack worked. Thank you. Now, I am getting the following future warning:

/usr/local/lib/python3.6/dist-packages/theano/tensor/ FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  result[diagonal_slice] = x

Do you know why this warning is being thrown and how I can remedy this?

Yeah it’s a classical and harmless warning. Usually we just get rid of it with:

import warnings
warnings.simplefilter(action="ignore", category=FutureWarning)


Can someone provide an example of what this would look like using tt.stack? I tried tt.concatenate(…) like this answer suggests but I’m getting a “TypeError: Join cannot handle arguments of dimension 0”. Since the OP replied using “tt.stack worked”, I’m wondering if @awareeye ran into the same problem?