How to set the priors for Bayesian estimation of Multivariate Normal Distribution when the correlation matrix has small values?

I am estimating the parameters for mean and covariance in Multivariate Normal Distribution (MVN). Following Wikipedia, I used MVN for the mean and Inverse-Wishart for covariance and tried Gibbs Sampling.

\begin{align} x_i &\sim \mathcal{N}(\mu,\Sigma), \ i=1,\ldots,N \end{align} Sampling $\mu$: \begin{align} p(\mu | x, \Sigma, \mu_0, \Sigma_0, \nu, \Psi) &\propto p(x | \mu, \Sigma) p(\mu | \mu_0, \Sigma_0) \ \cdots\ (1) \end{align} Sampling $\Sigma$: \begin{align} p(\Sigma | x, \Sigma, \mu_0, \Sigma_0, \nu, \Psi) &\propto p(x | \mu, \Sigma) p(\Sigma | \nu, \Psi) \ \cdots\ (2) \end{align}

Right hand side of (1) is MVN $\times$ MVN, so I draw a new $\mu$ in the $t$ th iteration following \begin{align} \mu^{(t)} &\sim \mathcal{N} (m, s)\\ s &= (\Sigma_{0}^{-1} + N \Sigma^{-1}_{(t)} )^{-1} \\ m &= s (\Sigma_{0}^{-1}\mu_0 + N \Sigma^{-1}_{(t)} \overline{x}) \\ \overline{x} &= \frac{1}{N} \sum_{i=1}^N x_i \end{align}

Right hand side of (2) is \begin{align} \Sigma^{(t)} &\sim InvW (N+\nu, \eta)\\ \eta &= \Psi + \sum_{i=1}^N (x_i - \mu^{(t)}) (x_i - \mu^{(t)})^T \end{align}

I implemented above in Python, but I could not recover the true values after enough number of iterations. Is this because of the priors?

JFYI, here is the Python code:

import numpy as np
import scipy as sp
import scipy.stats as sps
import numpy.random as npr
import seaborn as sns
import matplotlib.pyplot as plt
npr.seed(225)

# Create Data from Multivariate Normal
N = 1000 # number of data
D = 2 # dimensions
max_mean = 0.8
max_cov = 0.15

mean_vec = npr.normal(max_mean/2, 1, D)
cov_mat = npr.uniform(max_cov/2, max_cov, (D, D))
data = npr.multivariate_normal(mean_vec, cov_mat, N)

iter_num = 10000
show_num = 9500

# Initialization
# Prepare priors
## mean
mu_0 = np.repeat(0, D)
cov_0 = np.diag(np.repeat(0.5, D))

## cov
nu = D + 1
psi = np.identity(D)

mean_itr = npr.uniform(0, max_mean*2, D)
cov_itr = npr.uniform(0.01, max_cov*2, (D, D))

# Iteration
mean_chain = []
cov_chain = []
for i in range(iter_num):
# Update mean
cov0_inv = np.linalg.inv(cov_0)
cov_inv = np.linalg.inv(cov_itr)
cov_tmp = np.linalg.inv( cov0_inv  + N * cov_inv )
mean_tmp = cov_tmp.dot( cov0_inv.dot(mu_0) + N * np.dot(cov_inv,  data.mean(axis=0)) )
mean_itr = npr.multivariate_normal(mean_tmp, cov_tmp, 1)

mean_chain.append(mean_itr[0])

# Update cov
data_demean = data - mean_itr
scale_tmp = psi + (data_demean.transpose()).dot(data_demean).sum(axis=0)

cov_itr = sps.invwishart.rvs(N-1, scale_tmp)

cov_chain.append(cov_itr)

mean_chain = np.array(mean_chain)
cov_chain = np.array(cov_chain)

# Make FIgures
dim = 1
sns.distplot(mean_chain[show_num: , dim], hist=True, kde=False)
plt.plot([mean_vec[dim], mean_vec[dim]], [0, (iter_num-show_num)*0.2], linewidth=2, color='red')

index = (1,1)
sns.distplot(cov_chain[show_num:, index[0], index[1]], hist=True, kde=False)
plt.plot([cov_mat[index[0], index[1]], cov_mat[index[0], index[1]]], [0, (iter_num-show_num)*0.2], linewidth=2, color='red')


It seems means are fine, but it overestimates the covariance matrix when values are small (above: one of the means, below: one of the covariance, red lines are true values).

Updated: Using Normal-Inverse-Wishart

# Prepare priors
## mean
mu_0 = np.repeat(0, D)
k0 = 0.1
cov_0 = np.diag(np.repeat(0.5, D))

## cov
v0 = D + 1.5
psi = (v0 - D - 1) * np.identity(D) # I am not sure what is proper

# Initialization
mean_itr = npr.uniform(0, max_mean*2, D)
cov_itr = sps.invwishart.rvs(v0, psi)

# Iteration
mean_chain = []
cov_chain = []
for i in range(iter_num):
# Update mean
data_mean = data.mean(axis=0)
mean_tmp = (k0*mu_0 + N * data_mean ) / (k0 + N)
k = k0 + N
mean_itr = npr.multivariate_normal(mean_tmp, cov_itr/k, 1)

mean_chain.append(mean_itr[0])

# Update cov
data_demean = data - mean_itr
scale_tmp = psi + (data_demean.transpose()).dot(data_demean) + (k0*N)/(k0+N) *
(data_mean-mu_0).transpose().dot(data_mean-mu_0)
v = v0 + N

cov_itr = sps.invwishart.rvs(v, scale_tmp)

cov_chain.append(cov_itr)

mean_chain = np.array(mean_chain)
cov_chain = np.array(cov_chain)


Looking at the sampling of the $\Sigma$, I think you might have forgotten to replace the mean $\mu$ with the sampled $\mu(t)$ from your posterior distribution of the mean.

So $\eta$ should be determined as

$\eta=\Phi+\sum_n^N(x_n-\mu(t))(x_n-\mu(t))^T$

Gibbs sampling is an iterative procedure--when you sample from the posterior distribution of a single variable, you need to replace the other variables on which the sampled variable is dependent with their previously sampled value.

• Thanks for your correction! I updated the code, but I still overestimate the covariance matrix. – user2978524 Nov 9 '17 at 16:54
• When you resample the mean, you also need to replace the variance with the previously sampled value. Use $\Sigma(t)$ for sampling of $\mu(t)$ too. – Karlsson Yu Nov 9 '17 at 22:11
• Also check the derivation of $\eta$ again. I am surprised that there isnt an $1/(N-1)$ term in front of the $\sum_n^N$. – Karlsson Yu Nov 10 '17 at 0:01
• Thanks again for your corrections. For the first point, I edited the question, but the code already uses updated values in each step, so probably there is another reason why IW does not work well for small correlation matrix. For the second point, I am not sure whether or not we need $1/(N-1)$. I checked Wikipedia and this note, and neither of them puts $1/(N-1)$ in $\eta$. Also, the results look better without $1/(N-1)$ when I test it with my code. – user2978524 Nov 10 '17 at 3:04
• No worries. Happy coding – Karlsson Yu Nov 10 '17 at 11:17