5
$\begingroup$

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 enter image description here can be used to define the covariance of the Gaussian proposal distribution."

"There are two obvious approaches: (1) an independence proposal enter image description here or (2) a random walk proposal enter image description here, where enter image description here 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.

enter image description here

This is what I got from my code. I chose a proposal that is fixed at enter image description here with a variance matrix of enter image description here. 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?

enter image description here

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?

enter image description here

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 enter image description here 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

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ 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$. $\endgroup$
    – Xi'an
    May 22, 2021 at 5:45
  • $\begingroup$ You can use inline math here to format your formulae. $\endgroup$
    – phipsgabler
    May 22, 2021 at 14:53
  • 1
    $\begingroup$ 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. $\endgroup$
    – Ian
    May 22, 2021 at 17:29
  • 1
    $\begingroup$ (Or, as it apparently turned out, an actual bug in the implementation of the desired algorithm.) $\endgroup$
    – Ian
    May 22, 2021 at 17:54

2 Answers 2

Reset to default
2
$\begingroup$

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))
}

enter image description here

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)

Sampled data from the grid

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$
    – Xi'an
    May 22, 2021 at 5:31
6
$\begingroup$

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:

enter image description here 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:

enter image description here

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.