# What is the correct form of Metropolis Hasting step in scaled Inverse Wishart prior for covariance matrix?

I was going through the paper of O'Malley and Zaslavsky (2008) for the scaled inverse Wishart priors for a covariance matrix, in order to write an R-code for hierarchical Bayesian estimation of mixed logit model. They have decomposed the covariance matrix as $\Sigma = \Delta\Phi\Delta$; $\Delta$ is a diagonal matrix of $\delta_i$ and $\log(\delta_i) \sim \mathcal{N}(m,s)$; $\Phi$ is a correlation matrix, $\Phi \sim IW(nu,I)$.

They have used MCMC methods, combination of Gibbs sampling with Metropolis Hastings steps. In order to update the $\delta_i$ they used Metropolis Hasting step where the proposal distribution is the logarithm of a t distribution with 3 degrees of freedom. The details are given in the following picture:

I have defined the metropolis hasting step as But I'm not sure whether this is correct or not, as I'm getting non-positive definite covariance matrices. Should I use the proposed distribution in Metropolis Hasting or not?

I found this question when I was working on inverse-Wishart in this question. I wrote a piece of code in Python. To do MH, I think we only need to compare old/new likelihood and prior.

# Load Libraries
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
import collections
npr.seed(225)

# Create Data
N = 100 # number of data
D = 5 # dimensions
max_mean = 0.8
max_cov = 0.15

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

# Likelihood Functions
def mvn_loglik_single(x, mu, cov):
loglik = -0.5 * ( (x - mu).dot( np.linalg.inv(cov) ).dot((x-mu).transpose()) ) - 0.5 * np.log( np.linalg.det(2.0 * np.pi * cov ) )
return loglik
def mvn_loglik(x, mu, cov):
num = x.shape[0]
if num == 1:
return mvn_loglik_single(x, mu, cov)
else:
loglik = 0.0
for i in range(num):
loglik += mvn_loglik_single(data[i], mean_vec, cov_mat)

return loglik

# Initialization
iter_num = 2500
show_num = 1800

# 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)
delta = np.identity(D)

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

# store_delta
delta_proposed = np.zeros((D,D))

# Iteration
# Iteration
mean_chain = []
cov_chain = []
accept = 0
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 delta (update all at the same time)
for k in range(D):
delta_proposed[k,k] = delta[k,k] * np.exp( npr.normal(0, 0.02) )
cov_proposed = delta_proposed.dot(sigma_itr).dot(delta_proposed)

diflikelihood = mvn_loglik(data, mean_itr, cov_proposed) - mvn_loglik(data, mean_itr, cov_itr)
for k in range(D):
# consider prior
diflikelihood += sps.gamma.pdf(delta_proposed[k,k], a=2.0, scale=2.0) - sps.gamma.pdf(delta[k,k], a=2.0, scale=2.0)

r = min(0, diflikelihood) # np.log(1) = 0
u = np.log(npr.uniform(0,1))

if u < r:
accept  += 1
delta = delta_proposed
else:
pass

# Update cov
data_demean = data - mean_itr
scale_tmp = psi + np.linalg.inv(delta).dot( (data_demean.transpose()).dot(data_demean) ).dot(np.linalg.inv(delta))

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

cov_itr = delta.dot(sigma_itr).dot(delta)
cov_chain.append(cov_itr)

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

accept / iter_num

# Figures
dim = 0
sns.pointplot(x=np.arange(show_num, iter_num, 1), y=mean_chain[show_num: , dim])
plt.plot([0, iter_num-show_num], [mean_vec[dim], mean_vec[dim]], linewidth=2, color='red')

index = (2,2)
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')


Here are figures for covariance matrix. Not all parameters are well estimated probably because I moved all parameters in $\Delta$ at once to speed up, not one by one, which makes the results worse when we have more than seven or eight dimensions.

• Code (see meaning 7) is an uncountable noun (/a mass noun, like "rice"); so "some code" or "a piece of code" rather than "a code" (or perhaps you could say something like "a macro", "a function" or "a program"). Nov 9, 2017 at 20:46
• Could you explain a bit what are the graphs? They are the histograms of what? Jul 3, 2019 at 10:40
• @Anoldmaninthesea. It's a histogram of sampled values. Sep 24, 2019 at 21:55