3

我想从后验采样，其中 LambdaA 和 LambdaB 是 A 和 B 的指数率。此外，y 是 rv 的观察值。

后验由下式给出

在此处输入图像描述

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

数据：

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 值

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) 进行比较。有什么帮助可以使此代码正常工作吗？

4

1 回答 1

4

如果你真的认为随机游走意味着

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

您应该重新考虑并投入阅读马尔可夫链理论和马尔可夫链蒙特卡洛的基础：在每次迭代中，您将一个统一 U(0,1) 变量添加到当前值。因此，您总是建议增加当前值。你认为这会产生一个遍历马尔可夫链吗？

还有一个错误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 中产生对称的结果。

于 2017-11-30T09:40:08.687 回答