# Metropolis Hastings algorithm bivariate normals

I need some help implementing the (1) independence Gaussian proposal and (2) random walk Gaussian proposal to simulate from a mixture bivariate normal distribution.

"If we have a continuous state space, the Hessian H at a local mode can be used to define the covariance of the Gaussian proposal distribution."

"There are two obvious approaches: (1) an independence proposal or (2) a random walk proposal , where where D is the dimensionality of w, resulting in an acceptance rate of .234."

This issue is also discussed here: https://forum.dynare.org/t/finding-the-right-draws-from-proposals-in-a-metropolis-hastings/12812 "In Dynare the default jumping distribution is the Gaussian distribution centred on the previous state of the chain and with a covariance matrix given by the inverse of the Hessian matrix of the posterior kernel evaluated at the posterior mode."

However, I do not know what the Hessian matrix at the posterior kernel evaluated at the posterior mode means. Here is my attempt at implementing the M-H algorithm.

Part I: Independence proposal

This is what the mixture bivariate distribution should look like.

This is what I got from my code. I chose a proposal that is fixed at with a variance matrix of . This only seems to get two mixtures. Needless to say this does not follow the instructions that I should use the Hessian at a certain point. What Hessian at the local mode mean?

library(MASS)
mu=c(2,2)
Sigma=matrix(c(1, 1/2*1*sqrt(2), 1/2*1*sqrt(2), 2), nrow=2)
mu2=c(2,8)
Sigma2=matrix(c(1, -1/2*1*sqrt(2), -1/2*1*sqrt(2), 2), nrow=2)
mu3=c(6,4)
Sigma3=matrix(c(3, -1/2*sqrt(3)*sqrt(2), -1/2*sqrt(3)*sqrt(2), 2), nrow=2)
dat=data.frame(matrix(0, nrow=0, ncol=2))
for (i in 1:3000) {
u=runif(1)
if (u<1/3) {
dat[i,]=mvrnorm(1, mu, Sigma)
} else if (u<2/3) {
dat[i,]=mvrnorm(1, mu2, Sigma2)
} else {
dat[i,]=mvrnorm(1, mu3, Sigma3)
}
}
library(ggplot2)
ggplot(dat, aes(X1, X2)) + stat_density_2d(aes(fill=..level..), geom="polygon", color="white")

d=function(x) {
a=1/(2*pi)/sqrt(det(Sigma))*exp(-.5*t(x-mu)%*%solve(Sigma)%*%(x-mu))
b=1/(2*pi)/sqrt(det(Sigma2))*exp(-.5*t(x-mu2)%*%solve(Sigma2)%*%(x-mu2))
c=1/(2*pi)/sqrt(det(Sigma3))*exp(-.5*t(x-mu3)%*%solve(Sigma3)%*%(x-mu3))
return(1/3*(a+b+c))
}

#Independence

sig=matrix(c(4, 1/2*2*2, 1/2*2*2, 4), nrow=2)
y=data.frame(matrix(0, nrow=0, ncol=2))
x=c(4,4)
accepteds=0
for (i in 1:10000) {
print(i)
xp=mvrnorm(1, c(4,4), sig)
a=d(xp)/d(x)
r=min(1, a)
u=runif(1)
if (u<r) {
x=xp
accepteds=accepteds+1
}
y[i,]=x
}
accepteds/10000
ggplot(y, aes(X1, X2)) + stat_density_2d(aes(fill=..level..), geom="polygon", color="white")

Is this the correct Hessian, and if so what mean and variance would I use?

hessian=function(x, mu, Sigma) {
return(1/(2*pi)/sqrt(det(Sigma))*exp(-1/2* (t(x-mu) %*% solve(Sigma) %*% (x-mu))[1,1]) * (-1/2*2*solve(Sigma)))
}

Part II: Random walk proposal

I have decided to keep this in the same question. The random walk proposal does not work as the turns out to be not positive definite (not valid variance matrix). Any idea what I could be doing wrong?

dmvrnorm=function(x, mu, Sigma) {
return(1/(2*pi)/sqrt(solve(Sigma))*exp(-.5*(t(x-mu)%*%solve(Sigma)%*%(x-mu))[1,1]))
}

sig=matrix(c(4, 1/2*2*2, 1/2*2*2, 4), nrow=2)
y=data.frame(matrix(0, nrow=0, ncol=2))
x=c(4,4)
accepteds=0
for (i in 1:10000) {
print(i)
sig=solve(hessian(x, x, sig))*2.38^2/2
print(sig)
xp=mvrnorm(1, x, sig)
a=d(xp)*dmvrnorm(x, xp, sig)/d(x)/dmvrnorm(xp, x, sig)
r=min(1, a)
u=runif(1)
if (u<r) {
x=xp
accepteds=accepteds+1
}
y[i,]=x
}
accepteds/10000
ggplot(y, aes(X1, X2)) + stat_density_2d(aes(fill=..level..), geom="polygon", color="white")

Error: Error in mvrnorm(1, x, sig) : 'Sigma' is not positive definite

• The Hessian is the matrix of the second derivatives of the log-density. In the mixture example, the mode appears to be near $\mu_2$. May 22 at 5:45
• You can use inline math here to format your formulae. May 22 at 14:53
• In principle, as long as every point is reachable from every other point, MH eventually converges. So your bad behavior must be caused by either too few iterations or too many rejections or both.
– Ian
May 22 at 17:29
• (Or, as it apparently turned out, an actual bug in the implementation of the desired algorithm.)
– Ian
May 22 at 17:54

Since the independent Metropolis-Hastings algorithm is formally valid, the issue stands in an inadequate calibration of the proposal to reach the entire support of the target (mixture) distribution. I just modified the code by choosing a larger variance matrix

sig=5*matrix(c(4, 1/2*2*2, 1/2*2*2, 4), nrow=2)

ran the chain 10⁵ iterations, and somewhat recovered the entire target:

However, there is an error in the code and possibly a misunderstanding of the independent Metropolis-Hastings algorithm. The acceptance probability

a=d(xp)/d(x)
r=min(1, a)

should divide the ratio of the targets by the ratio of the proposals

a=d(xp)/d(x)/dmvnorm(xp,c(4,4), sig)*dmvnorm(x,c(4,4), sig)
r=min(1, a) #superfluous for the acceptance

which returns a better representation of the target after 10⁴ iterations:

EDIT: REFER TO THE CORRECT ANSWER ABOVE:

Note that this is a problem of MIXTURE of NORMALS rather than BIVARIATE NORMALS:

The issue here is that the proposal density does not cover the whole region of interest. You were to calibrate it accordingly to cover the ROI as explained in the answer above. Here is another example of using a grid/uniform distribution ie

we can see from the graph that the x and y axis runs from -2 to 12. Then you could use that to propose. and since it is symmetric/constant, we can simply implement metropolis rather than metropolis Hastings.

library(mvtnorm)
mu <- c(2,2)
Sigma <- matrix(c(1, 1/2*1*sqrt(2), 1/2*1*sqrt(2), 2), nrow=2)
mu2 <- c(2,8)
Sigma2 <- matrix(c(1, -1/2*1*sqrt(2), -1/2*1*sqrt(2), 2), nrow=2)
mu3 <- c(6,4)
Sigma3 <- matrix(c(3, -1/2*sqrt(3)*sqrt(2), -1/2*sqrt(3)*sqrt(2), 2), nrow=2)

n <- 3000
s <- rmultinom(1, n, c(1,1,1))
dat <- mapply(rmvnorm, s, mean = list(mu, mu2, mu3), sigma = list(Sigma, Sigma2, Sigma3))
dat1 <- setNames(do.call(rbind.data.frame, dat), c("X1", "X2"))

ggplot(dat1, aes(X1, X2)) +
stat_density_2d(aes(fill=..level..), geom="polygon", color="white")

d <- function(x) {
a <- dmvnorm(x, mu, Sigma)
b <- dmvnorm(x, mu2, Sigma2)
c <- dmvnorm(x, mu3, Sigma3)
mean(c(a,b,c))
}

## Implementation of Metropolis:

B <- 40000
y <- data.frame(matrix(nrow = B, ncol = 2))
colnames(y) <- c("X1", "X2")
y[1, ] <-colMeans(dat1)
accept<- 1
for(i in seq(2,B)){
prop <- c(runif(1, -2, 12),runif(1, -2,12))#
if( runif(1)<d(prop)/d(y[i-1, ])) {
y[i, ]<- prop
accept <- accept + 1
}
else y[i, ]<- y[i-1, ]
}

##Burn the first 5000 points
ggplot(y[-seq(5000), ], aes(X1, X2)) +
stat_density_2d(aes(fill=..level..), geom="polygon", color="white")
print(accept/B)

• I disagree with the statement that one "cannot propose from a normal density to capture a mixture of normals" since the Metropolis-Hastings algorithm attached to a normal proposal returns an ergodic Markov chain with the proper stationary distribution. May 22 at 5:31