我需要一些帮助来实现（1）独立高斯建议和（2）随机游走高斯建议，以从混合二元正态分布中进行模拟。

“如果我们有一个连续的状态空间，则可以使用局部模式下的 Hessian H来定义高斯提议分布的协方差。”

“有两种明显的方法：（1）独立提议或（2）随机游走提议，其中D是w的维数，导致接受率为 0.234。”

这里也讨论了这个问题：https : //forum.dynare.org/t/finding-the-right-draws-from-proposals-in-a-metropolis-hastings/12812 “在 Dynare 中，默认的跳跃分布是高斯分布分布以链的前一个状态为中心，协方差矩阵由在后验模式下评估的后验核的 Hessian 矩阵的逆矩阵给出。”

但是，我不知道在后验模式下评估的后验核处的 Hessian 矩阵是什么意思。这是我实现 MH 算法的尝试。

第一部分：独立提案

这就是混合双变量分布的样子。

这是我从我的代码中得到的。我选择了一个固定为方差矩阵的提案。这似乎只得到两种混合物。不用说，这不符合我应该在某个时候使用 Hessian 的说明。本地模式下的 Hessian 是什么意思？

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

这是正确的 Hessian，如果是，我将使用什么均值和方差？

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

第二部分：随机游走提案

我决定将其保留在同一个问题中。随机游走提议不起作用，因为结果不是正定的（不是有效的方差矩阵）。知道我做错了什么吗？

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 in mvrnorm(1, x, sig) : 'Sigma' is not positive definite