R - 使用 gibbs 的 RW metropolis 失败

Question

我想从后验中抽样，其中 LambdaA 和 LambdaB 是 A 和 B 的指数率。此外，y 是 r.v.'s.

的观测值

后验由

给出

出于数字原因，我正在记录此函数。

数据：

n<-100
y<- c(rexp(n))

后验对数：

dmix<-function(LambdaA,LambdaB,w){


ifelse( LambdaA<=0|LambdaB<=0|w<0|w>1 ,0,log(w*LambdaA*LambdaB*exp(-2*(LambdaA+LambdaB))*prod(w*LambdaA*exp(-
LambdaA*y) + (1-w)*LambdaB*exp(-LambdaB*y)) ))}

U-values

U.lambdaB <- runif(1)
U.lambdaA<- runif(1)
U.w<- runif(1)

数步数

REJLambdaB <- 1
REJw <- 1
REJLambdaA<-1

初始点数

LambdaB <- LambdaA<- w<- numeric(n)
LambdaA[1]<-0.5
LambdaB[1] <- 0.5
w[1] <- 0.5

随机游走MH算法，一次更新每个分量：

for (t in 2:n){


 LambdaBprop<- rnorm(1,LB[t-1],0.5)
 wprop<- rnorm(1,w[t-1],0.5)
 LambdaAprop<- rnorm(1,LB[t-1],0.5)

logalpha1 = dmix(LambdaAprop,LambdaB[t-1],w[t-1])-dmix(LambdaA[t-1],LambdaB[t-
1],w[t-1])

logalpha2 = dmix(LambdaA[t-1],LambdaBprop,w[t-1])-dmix(LA[t-1],LB[t-1],w[t-
 1])




if (!is.null(log(U.lambdaB) > logalpha2))

{LambdaB[t] <- LambdaBprop}  ## accepted 

else{LambdaB[t] <- LambdaB[t-1] ##rejected
REJLambdaB<-REJLambdaB+1} 


if (!is.null(log(U.lambdaA) > logalpha1)) 

{LambdaA[t]<-LambdaAprop}

else {LambdaA[t]<-LambdaA[t-1]
REJLambdaA<-REJLambdaA+1}

 if (w[t]<0|w[t]>1) 

{w[t]<-w[t-1]}

else {w[t]<-wprop
REJw<-REJw+1}

}

最终，我的后验出现问题，因为在评估 logalpha 时我总是得到无穷大或 0。请注意，我正在寻找比较 log($\alpha(x'|x))$ 和 log(U)。任何帮助使此代码工作？

Answer 1

如果你真的认为 random walk 意味着

lambdB[t]<- lambdB[t-1] + runif(1)
w[t]<- w[t-1] + runif(1)
lambdA[t] <- lambdB[t-1] + runif(1)

您应该重新考虑并投入阅读马尔可夫链理论和马尔可夫链的基础知识 Monte Carlo：在每次迭代中，您向当前值添加一个统一的 U(0,1) 变量。因此，您总是建议增加当前值。您认为这会产生 ergodic Markov chain 吗？

dmix中还有一个错误：因为你使用的是对数，所以记住log(0)=-oo。数量 logalpha1 和 logalpha2 未正确更新。还有更多的编程错误，比如 !is.null 的不正确使用...无论如何，这是一个有效的更正 R 代码：

n<-100
y<- c(rexp(n))

#Logarithm of posterior:

dmix<-function(LambdaA,LambdaB,w){

  ifelse( (LambdaA<=0)|(LambdaB<=0)|(w<0)|(w>1) ,
    -1e50,log(w*LambdaA*LambdaB)-2*(LambdaA+LambdaB)+sum(log(w*LambdaA*exp(-
    LambdaA*y) + (1-w)*LambdaB*exp(-LambdaB*y))) )}

#Count steps

REJLambdaB <- 1
REJw <- 1
REJLambdaA<-1

#Initial points
N <- 1e4
LambdaB <- LambdaA <- w<- numeric(N)
LambdaA[1] <- LambdaB[1] <- w[1] <- 0.5

U.lambdaB <- runif(N)
U.lambdaA<- runif(N)
U.w <- runif(N)

for (t in 2:N){
LambdaBprop=rnorm(1,LambdaB[t-1],0.5)
LambdaAprop=rnorm(1,LambdaA[t-1],0.5)
wprop=rnorm(1,w[t-1],0.05)

logalpha2 = dmix(LambdaA[t-1],LambdaBprop,w[t-1])-dmix(LambdaA[t-1],LambdaB[t-1],w[t-1])

if ((log(U.lambdaB[t]) < logalpha2))
  {LambdaB[t] <- LambdaBprop}  ## accepted
else{LambdaB[t] <- LambdaB[t-1] ##rejected
     REJLambdaB<-REJLambdaB+1}

logalpha1 = dmix(LambdaAprop,LambdaB[t],w[t-1])-dmix(LambdaA[t-1],LambdaB[t],w[t-1])

if ((log(U.lambdaA[t]) < logalpha1)) 
  {LambdaA[t]<-LambdaAprop}
else {LambdaA[t]<-LambdaA[t-1]
  REJLambdaA<-REJLambdaA+1}

logw = dmix(LambdaA[t],LambdaB[t],wprop)-dmix(LambdaA[t],LambdaB[t],w[t-1])

if (w[t]<0|w[t]>1|(log(U.w[t])>logw))
    {w[t]<-w[t-1]}
else {w[t]<-wprop
    REJw<-REJw+1}

}

如结果所示

后验在 Lambda 中产生对称结果。

R - 使用 gibbs 的 RW metropolis 失败

R - RW metropolis using gibbs fails

mcmc

random-walk