# Migration from pymc2 with deterministic function on matrixes

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

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
``````

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**2, rho*sigmas*sigmas], [rho*sigmas*sigmas, sigmas**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

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

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));
`````` 2 Likes

Thanks again!