R中的这两个随机游走Metropolis-Hastings算法实现之间有什么区别吗？
Is there any difference between these two random walk Metropolis-Hastings algorithm implementations in R?
我正在尝试在 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 将所有内容分解为单独的部分，并按顺序进行。
显然，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
问题是您使用的是随机数，因此要获得可重现的结果，您需要在 运行 算法之前使用 set.seed。
我从 for 循环中取出函数的定义并使用 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
我正在尝试在 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 将所有内容分解为单独的部分，并按顺序进行。
显然，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
问题是您使用的是随机数，因此要获得可重现的结果，您需要在 运行 算法之前使用 set.seed。
我从 for 循环中取出函数的定义并使用 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