I have the following code,
gibbs <-function(m,theta = 0.25, lambda =0.55, n =1){
    alpha <- 1.5
    beta <- 1.5
    gamma <- 1.5
    x<- array(0,c(m+1, 3))
    x[1,1] <- theta
    x[1,2] <- lambda
    x[1,3]<- n
    for(t in 2:(m+1)){
        x[t,1] <- rbinom(1, x[t-1,3], x[t-1,1])
        x[t,2]<-rbeta(1, x[t-1,1] + alpha, x[t-1,3] - x[t-1,1] + beta)
        x[t,3] <- rpois(1,(1 - x[t-1,1])*gamma)
    }
    x
}
gibbs(100)

it returns only 1 or 2 values of theta, some times NaN, this may be if any
theta is greater than 1, which is used as the probability for the next
rbinom(), so returns NaN. Can some one suggest to solve this problem?

	[[alternative HTML version deleted]]

