[R] Gibbs sampler

[Duncan Murdoch]

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.