I’m trying to implement what seems like a fairly simple model, and am having some issues.

I’m pretty new to pymc3 – any help would be really appreciated!

Cheers

David

### Data

Split into groups, indexed by i:

- a_i: shape (n_i, 2)
- b_i: shape (m_i, 2)
- y_i

### Model

Linear model with:

y_i \sim N(\alpha + \beta x_i, \sigma)

Where x_i comes from the euclidean distance between samples of two multivariate normal distributions:

x_i = ||a_i - b_i||

a_i \sim MN(\mu_{ai}, cov_{ai})

b_i \sim MN(\mu_{bi}, cov_{bi})

Using uninformative priors.

### Issue

I *think* what needs to be programmed, is that for every sampling step:

- values are drawn from the two multivariate normal distributions, which have a_i and b_i as observed data, respectively.
- x_i is computed from the values that are drawn.
- N(y_i | \alpha + \beta x_i, \sigma) is computed.

What I have tried to implement won’t do this. For one thing the function `mvNormal`

does not return a sample from the instance of `pm.MvNormal`

. But these `pm.MvNormal`

instances do not have a `.random()`

method because they have observed data.

It seems like this would be possible if \mu_{ai}, cov_{ai}, \mu_{bi}, and cov_{bi} were estimated in a separate step, but that there should be a way of combining it all into a single model.

Does this model make sense? And if so, how should it be implemented?

### Attempted implementation

```
import pymc3 as pm
import numpy as np
def mvNormal(observed, pre):
"""Model observed data with a multivariate normal distribution.
Args:
data (np.array): Shape [n, 2].
pre (str): Prefix for variable names.
Returns:
pm.MvNormal
"""
mu = pm.Normal(pre + 'mu', mu=np.zeros(2), sd=np.ones(2) * 10, shape=(2,))
sd_dist = pm.HalfCauchy.dist(beta=2.5, shape=(2,))
chol_packed = pm.LKJCholeskyCov(pre + 'cp', n=2, eta=2, sd_dist=sd_dist)
chol = pm.expand_packed_triangular(2, chol_packed)
return pm.MvNormal(pre + 'vals', mu=mu, chol=chol, observed=observed)
def distsMvNormals(a, b, pre):
"""Distribution of distances between samples from 2 multivariate normal
distributions.
Args:
a (np.array): Observed data. Shape [n, 2].
b (np.array): Observed data. Shape [m, 2].
pre (str): Prefix for variable names.
Returns:
Euclidean distance between samples of a and b.
"""
var_a = mvNormal(a, pre + 'mvn-a-')
var_b = mvNormal(b, pre + 'mvn-b-')
return pm.Deterministic(pre + 'dist', (var_a - var_b).norm(L=2))
with pm.Model() as model:
n, m = 50, 50
k = 3
x = [
distsMvNormals(np.random.randn(n, 2), np.random.randn(m, 2), "a-"),
distsMvNormals(np.random.randn(n, 2), np.random.randn(m, 2), "b-"),
distsMvNormals(np.random.randn(n, 2), np.random.randn(m, 2), "c-")
]
alpha = pm.Normal('alpha', mu=0, sd=10)
beta = pm.Normal('beta', mu=0, sd=10)
sigma = pm.HalfNormal('sigma', sd=10)
y = pm.Normal('y', mu=alpha + beta * x, sd=sigma, observed=np.random.randn(k))
pm.sample(100) # Low n. draws for testing
```