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? 
 A: 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.
 
