ギブスを使用したR - RWメトロポリスは失敗します

Question

私は後方からサンプリングしたいと思います。ここでLambdaAとLambdaBはAとBの指数レートです。また、yはr.v.の観測値です。

後部は

によって与えられ、数値的な理由から、私はこの機能のログを取っています。

データ：後方の

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) 

Countが

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

ステップ

初期点

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

ランダムウォークMHアルゴリズム、一度に各コンポーネントを更新：私はlogalpha年代を評価する際に無限大または0のいずれかを取得しておくので、

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} 

} 

は最終的に、私は私の後部に問題が生じています。ログ log（$ \ alpha（x '| x））$とlog（U）を比較することを検討しています。このコードを動作させる助けとなることはありますか？あなたが本当にrandom walkが

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

を意味すると思われる場合

Answer 1

あなたは再考し、マルコフ連鎖の理論とマルコフ連鎖モンテカルロの拠点を読みに投資する必要があります。各反復で、あなたは制服U（0,1を追加します）を現在の値に変更します。したがって、あなたは常にに提案するを増やす現在の値。これでergodic Markov chainが生成されると思いますか？

dmixにも間違いがあります：対数を扱うので、log（0）= - ooを覚えておいてください。また、数量logalpha1とlogalpha2は正しく更新されません。

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} 

} 

は後部が生成する結果が示すように：と、より多くのプログラミングエラーは、!is.nullの誤った使用のような...とにかく、ここで働く修正Rコードがありますラムダの対称的な結果。