R中这两个随机游走Metropolis-Hastings算法实现之间有什么区别吗？ [关闭]

Question

我正在尝试在 R 中实现随机游走 Metropolis-Hastings 算法。我使用了自定义函数 logit 和 invlogit 来应用和撤消 logit 函数。我还使用正态分布来添加随机噪声。鉴于这两个事实，一旦你对变换后的参数 + 随机噪声使用逆变换，建议分布不再是对称的，这就是为什么我随后将校正项 log(yt*(1 - yt)) - log(xt*(1 - xt)) 应用于接受概率。

我的问题是，似乎有两种方法可以在 R 中实现此算法。如果这两种方法相等，那么，据我所知，在计算 acceptanceRate 时，我应该得到相等的值。然而，事实并非如此，这让我相信一种实现有缺陷（有错误），而另一种则没有。

但是，另外两种可能性是（1）两种方式都不正确或（2）两种方式都是正确的，我误解了一些东西。我是 R 编码的新手，我仍然无法理解为什么两个实现之间的 acceptanceRate 的值存在这些偏差。

注意：我感兴趣的具体问题是为什么我在两个实现之间得到不同的 acceptanceRate 值。

实施 1

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)

B <- 10000           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point
  logit <- function(p) log(p/(1-p))
  invlogit <- function(lo) 1/(1 + exp(-lo))
  yt <- invlogit(rnorm(1, mean = logit(xt), sd = 0.45))      ## proposal
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate <- num.accept/B

看看实现 1 如何使用yt <- invlogit(rnorm(1, mean = logit(xt), sd = 0.45))？一切都是积累和共同完成的。

实施 2

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)

B <- 10000           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point
  logit <- function(p) log(p/(1-p))
  xt <- logit(xt)
  yt <- xt + rnorm(1, mean = 0, sd = 0.45)      ## proposal
  invlogit <- function(lo) 1/(1 + exp(-lo))
  xt <- invlogit(xt)
  yt <- invlogit(yt)
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate <- num.accept/B

请注意，实现 2 将所有内容分解为单独的部分，按顺序进行。

Answer 1

问题在于您使用的是随机数，因此要获得可重现的结果，您需要在运行算法之前使用set.seed。
我从for-loop 中取出函数的定义并使用set.seed。在这两种情况下我得到了相同的结果：

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)
logit <- function(p) log(p/(1-p))
invlogit <- function(lo) 1/(1 + exp(-lo))

第一次实施

set.seed(42)
B <- 10000  ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point
  yt <- invlogit(rnorm(1, mean = logit(xt), sd = 0.45))      ## proposal
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate1 <- num.accept/B

rm(B, chain, num.accept, i, lapt, xt, yt)

第二次实施

set.seed(42)
B <- 10000           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals

for(i in 1:B){
  xt <- chain[i] ## current point
  xt <- logit(xt)
  yt <- xt + rnorm(1, mean = 0, sd = 0.45)      ## proposal
  xt <- invlogit(xt)
  yt <- invlogit(yt)
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate2 <- num.accept/B

acceptanceRate1
# [1] 0.7029
acceptanceRate2
# [1] 0.7029

Answer 2

显然，OP 在没有设置种子的情况下比较了依赖随机数生成器的两个函数 (set.seed)。

我看不出它有什么问题。对于小型连锁店，我得到相同的结果。

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)
invlogit <- function(lo) 1/(1 + exp(-lo))
logit <- function(p) log(p/(1-p))

set.seed(1)

B <- 100           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point

  xt <- logit(xt)
  yt <- xt + rnorm(1, mean = 0, sd = 0.45)      ## proposal

  xt <- invlogit(xt)
  yt <- invlogit(yt)
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate <- num.accept/B

# acceptanceRate 
# [1] 0.69

# chain[30:40]  
# [1] 0.7674114 0.6612332 0.5867199 0.5867199 0.5744098 0.6033942 0.5359917  [8] 0.5359917 0.5359917 0.6040635 0.6040635

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)
logit <- function(p) log(p/(1-p))
invlogit <- function(lo) 1/(1 + exp(-lo))

set.seed(1)
B <- 100           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point

  yt <- invlogit(rnorm(1, mean = logit(xt), sd = 0.45))      ## proposal
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate <- num.accept/B

# acceptanceRate 
# [1] 0.69

# chain[30:40]  
# [1] 0.7674114 0.6612332 0.5867199 0.5867199 0.5744098 0.6033942 0.5359917  [8] 0.5359917 0.5359917 0.6040635 0.6040635