R help - Nov 2011 - Gibbs sampler

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]]

On 11-11-10 1:28 PM, Gyanendra Pokharel wrote:
> 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)
Are you certain that x[t-1,3] - x[t-1,1] + beta will always be positive? 
  If it is negative, rbeta will return NaN.

Duncan Murdoch
>          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]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.