Migration from pymc2 with deterministic function on matrixes

Questions
1

Apologies in advance for the imprecise language. I am attempting to sample a covariance matrix from a dataset with two variables X and Y (both standardized). When migrating from pymc2 to pymc3 the definition of the prior and models is easy:

with pm.Model() as model:
    # assume an uniformative prior for rho between -1 and 1 and starting at 0 for the sampling
    rho = pm.Uniform('rho', lower=-1, upper=1)
    # the priors on our expected data values and starting point are derived from the data
    mu1 = pm.Normal('mu1', 0.0, 1.0/(1000)**2.0, observed = x.mean())
    mu2 = pm.Normal('mu2', 0.0, 1.0/(1000)**2.0, observed = y.mean())
    sigma1 = pm.InverseGamma('sigma1', 11.0, 10.0, observed = x.std())
    sigma2 = pm.InverseGamma('sigma2', 11.0, 10.0, observed = y.std())

but while in pymc2 I could define two deterministic functions

@pm.deterministic
def mean(mu1=mu1, mu2=mu2):
    return np.array([mu1, mu2])    
    
@pm.deterministic
def CoV(rho=rho, sigma1=sigma1, sigma2=sigma2):
    ss1, ss2, rss = sigma1 * sigma1, sigma2 * sigma2, rho * sigma1 * sigma2
    return [[ss1, rss], [rss, ss2]]

xy = pm.MvNormalCov('xy', mu=mean, C=CoV, value=data.T, observed=True)

Now the equivalent

 mean_val = pm.Deterministic('mean value', [mu1, mu2])
 CoV = pm.Deterministic('Covariance', [[sigma1**2.0, rho*sigma1*sigma2],[rho*sigma1*sigma2, sigma2**2.0]])
 xy = pm.MvNormal('xy', mu=mean_val, cov=CoV, observed=data.T)

triggers an error

 mean_val = pm.Deterministic('mean value', np.array([mu1, mu2]))
  File "/scratch/home/pfigueir/anaconda3/envs/py37/lib/python3.7/site-packages/pymc3/model.py", line 1597, in Deterministic
    var = var.copy(model.name_for(name))
TypeError: order not understood

That I think I should address by reformulating the question. Could you please help me?

2

Hi pfigueira

Here is code that can be used in PyMC3.

N = 256
true_mus = np.array([-1, 3])
true_cov = np.array([[1.4*1.4, 1.4*0.7*-0.2], [1.4*0.7*-0.2, 0.7*0.7]])
data = np.random.multivariate_normal(true_mus, true_cov, size=N)

with pm.Model() as model:
    mu = pm.Normal('mu', 0., 5., shape=2)
    rho = pm.Uniform('rho', lower=-1., upper=1.)
    sigmas = pm.InverseGamma('sigmas', 11., 10., shape=2)
    cov = tt.stack(([sigmas[0]**2, rho*sigmas[0]*sigmas[1]], [rho*sigmas[0]*sigmas[1], sigmas[1]**2]))
    obs = pm.MvNormal('obs', mu=mu, cov=cov, observed=data)

with model:
    trace = pm.sample()
print(pm.summary(trace))

In general though I would recommend using the LKJ Cholesky distribution when working with covariance matricies, especially when number of dimensions is above 2.

Another point, you do not need to specify the observed attributes for each of the mus and sigmas. Just passing in observed data on the MvNormal line will suffice.

2 Likes
3

Hello nkaimcaudle,

Thanks a lot, that was extremely helpful!

There are two little questions that remain:

  • You generate data by drawing N times from the true_mus and true cov, that were derived from (X,Y). Is there a way to feed directly the (X, Y) points to the MvNormal distribution?
    *I read the information you sent about LKJ Cholesky distribution with much interest. I just find it a bit hard to visualize how to include my uniform prior on rho. From this post I would say that it is with eta=1. Any comments?

Thanks a lot for your help!

Pedro

4

I don’t quite understand your question. Are (X, Y) mu parameters? If so then yes you can fix the mu parameters in the model easily.

mu = tt.stack((-0.9, 2.95))

eta=1 would indeed be a uniform prior on rho (-1, 1)

This code will let you visualise the effect of eta prior on rho

with pm.Model() as m_slope:
    # cov matrix:
    eta = 1.
    sd_dist = pm.Exponential.dist(1)
    packed_chol = pm.LKJCholeskyCov('chol_cov', eta=eta, n=2, sd_dist=sd_dist)
    chol = pm.expand_packed_triangular(2, packed_chol, lower=True)
    
    # Extract rho and sds:
    cov = pm.math.dot(chol, chol.T)
    sigma_ab = pm.Deterministic('sigma_district', tt.sqrt(tt.diag(cov)))
    corr = tt.diag(sigma_ab**-1).dot(cov.dot(tt.diag(sigma_ab**-1)))
    r = pm.Deterministic('Rho', corr[np.triu_indices(2, k=1)])
    
    prior = pm.sample_prior_predictive(samples=10000)
plt.hist(prior['Rho'], bins=20, density=True);
plt.title('eta = {}'.format(eta));

2020-05-20_08-21

2 Likes
5

Thanks again!