Metropolis hastrings estimation of covariance matrix

Metropolis_hastings.R

# Title     : Multivariate Metropolis Hastings for covariance matrix estimation
# Objective : Estimate distribution of multivariate covariances
# Created by: jonlachmann
# Created on: 2020-05-17

# Jeffreys prior, log transformed
jeffreys_log = function (sigma) {
	-(nrow(sigma)+1)/2 * log(det(sigma))
}

# Likelihood for multivariate normal, log transformed
likelihood_multi_log = function (data, sigma) {
	sigma_inv = solve(sigma)
	mu = colMeans(data)
	exponent = -1/2 * sum(apply(data, 1, function(row) t(row-mu)%*%sigma_inv%*%(row-mu)))
	-nrow(data)/2 * log(det(sigma)) + exponent
}

# Posterior, log transformed
posterior_multi_log = function (data, sigma) {
	jeffreys_log(sigma) + likelihood_multi_log(data, sigma)
}

# Proposal generator for covariance parameters (cov matrix in, cov matrix out)
proposal_rand_multi = function (current) {
	p = nrow(current)
	for(y in 1:(p-1)) {
		for (x in (y+1):p) {
			current[y,x] = current[x,y] = runif(1,current[x,y]-0.07,current[x,y]+0.07)
		}
	}
	return(current)
}

# Metropolis Hastings algorithm, multivariate
metropolis_multi = function (start, proposal, density, data, iter=1000) {
  	samplevec = vector(mode="list", length=iter)
    samplevec[[1]] = start
	for (i in 2:iter) {
		propose = proposal(samplevec[[i-1]])
		ratio = exp(min(0,(density(data, propose) - density(data, samplevec[[i-1]]))))
	  	if (runif(1) < ratio) samplevec[[i]] = propose
	  	else samplevec[[i]] = samplevec[[i-1]]
        if (i %% 100 == 0) print(i)
  	}
  	return(samplevec)
}

# Generate data to use
N=1000
data3=rmvnorm(N, c(0,0,0), matrix(c(1,0.2,0.4,0.2,1,0.6,0.4,0.6,1),3,3))

# Starting values for covariance matrix
sigma_null = matrix(c(1,0,0, 0,1,0, 0,0,1),3,3)

# Generate the sample
sample_mh = metropolis_multi(sigma_null, proposal_rand_multi, posterior_multi_log, data3, iter=1000)

# Plot the results
dev.off()
par(mfrow=c(3,1))
sample_mh = unlist(sample_mh)
rho1 = sample_mh[seq(2, length(sample_mh), 9)]
rho2 = sample_mh[seq(3, length(sample_mh), 9)]
rho3 = sample_mh[seq(6, length(sample_mh), 9)]

plot(rho1, type="l")
plot(rho2, type="l")
plot(rho3, type="l")