R - RW metropolis using gibbs fails

Question

I want to sample from the posterior, where LambdaA and LambdaB are exponential rates of A and B. Also, y is the observations of the r.v.'s.

The posterior is given by

and for numerical reasons, i am taking the log of this function.

Data:

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

Logarithm of posterior:

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)

Count steps

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

Initial points

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

Random walk MH algorithm, updating each component at a time:

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}

}

Ultimately, I am having problems with my posterior since I keep getting either infinity or 0's when evaluating logalpha's. Note that i am looking to compare log($\alpha(x'|x))$ with log(U). Any help to get this code to work ?

Answer 1

If you really think that a random walk means

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

you should reconsider and invest into reading the bases of Markov chain theory and Markov chain Monte Carlo: At each iteration you add a Uniform U(0,1) variate to the current value. Therefore you always propose to increase the current value. Do you think this could ever produce an ergodic Markov chain?

There is also a mistake in dmix: since you work with the logarithm, remember that log(0)=-oo. And the quantities logalpha1 and logalpha2 are not updated correctly. And many more programming errors, like the incorrect use of !is.null... Anyway here is a corrected R code that works:

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}

}

As shown by the outcome

the posterior produces a symmetric outcome in the Lambda's.