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

A: 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.
Your Code:
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)


