Multivariate Normal - cannot recover input values in end-to-end simulation

1

I cannot recover the correct parameters from an end-to-end simulation. I am assuming I have data distributed 5D-Gaussian like this:

mus = np.array([4,2,3,1,6]) #mean values
K = mus.shape[0]
N = 100 # Total amount of data
sigmas = np.eye(5)+0.3*np.random.random((5,5)) #covariance matrix, not positively defined but OK
data = np.random.multivariate_normal(mus, sigmas, size=N) # data generation

I propose the same model for inference

with pm.Model() as modelCorr:

    # Prior over component means
    mus = pm.MvNormal('mu_', mu=pm.floatX(np.zeros(K)), tau=pm.floatX(0.1 * np.eye(K)), shape=(K,))

    # Cholesky decomposed LKJ prior over component covariance matrices
    L = pm.LKJCholeskyCov('packed_L_', n=K, eta=2., sd_dist=pm.HalfCauchy.dist(1))[0]

    # Convert L to sigma for convenience
    sigma = pm.Deterministic('sigma_' ,L.dot(L.T))

    # Specify the likelihood
    obs = pm.MvNormal('like', mu=mus, chol=L, observed=data)

The model compiles well and I can sample using

with modelCorr:
    traceCorr = pm.sample(chains=8)

Then I estimate the PPCs using,

with modelCorr:
    est = pm.sample_posterior_predictive(traceCorr, var_names=['mu_', 'packed_L_'])

and finally I compute the expected distribution of the data using,

x_mu = np.asarray(est.posterior_predictive.mu_)
x_L = np.asarray(est.posterior_predictive.packed_L_)

# reshape to stack all the chains
xx_mu = x_mu.reshape(-1, x_mu.shape[-1])
xx_L = x_L.reshape(-1, x_L.shape[-1])

# data estimation
N = xx_mu.shape[0]
estimation = np.zeros(xx_mu.shape)
for n in range(N):
    estimation[n, :] = pm.MvNormal.dist(mu=xx_mu[n,:], chol=pm.expand_packed_triangular(5, xx_L[n,:], lower=True)).eval()

The mean values of the posterior means are already very bad (array([-0.06958857, -0.01141641, -0.00285485, -0.00243548, 0.01582915])). And of course, the estimated values have a really wide spread.

2

Can you try and see if

az.plot_posterior(traceCorr, var_names=["mu_"])
az.plot_posterior(traceCorr, var_names=["sigma_"])

gives the expected posteriors? Also why don’t you simulate your sigma to be positive definite by

sigma = sigma*sigma^T

1 Like
3

It works indeed! The pdfs and mean values of the traceplots are correct, and if I use

> x_mu = np.asarray(traceCorr.posterior.mu_)
> x_L = np.asarray(traceCorr.posterior.packed_L_)

Mean values are now OK. So the problem arises when I synthetize the observed values from the new samples from the PPC. When comparing histograms from est.observed_data.like they are similar to the real data, so ppc is working. Maybe there is a problem with pm.expand_packed_triangular(5, xx_L[n,:], lower=True)?

4

So the normal use for sample_posterior_predictive is to generate the posterior predictive distribution which is marginalizing your model p(X | param) over the param variable using your posterior p(param | data) as your pdf for parameter values (prior predictive check would be similar but you weigh your parameters by your priors p(param)). As you observed, if you want to get samples from the posterior distribution directly you can use

traceCorr.posterior["mu_"]

That being said, I do not know what happens internally when you supply var_names=[“_mu”] to posterior_predictive_check. When you look at the distribution of _mu from here, it looks like it is just sampling from the prior. On the other hand the following

with modelCorr:
    sampled_mu = pm.Deterministic("mu_samples", mus)
    est = pm.sample_posterior_predictive(traceCorr, var_names=['mu_samples'])

seems to produce mu_samples which have the expected posteriors. So maybe a dev could tell what is going on here internally when you supply other parameters as var_names to sample_posterior_predictive.

5

The behavior of sample_posterior_predictive with var_names is a bit tricky. We have updated the documentation (not yet released) in Add examples explaining advanced applications of `sample_posterior_predictive` by ricardoV94 · Pull Request #7014 · pymc-devs/pymc · GitHub

We want to make the behavior more intuitive: Make `sample_posterior_predictive` API less surprising · Issue #7069 · pymc-devs/pymc · GitHub

In your example you probably don’t want to include it in var_names

6

Blockquote

Very interesting reading. As I understand ppc now, the only valid case for me will be to sample from observed variables (to obtain out of sample predictions).

Thank you both for the feedback!

3 Likes