model;
{   
   for (i in 1:N){
      p[i] <- 1/N
   }
   S ~ dcat(p[])
   sprime <- Rprime-rprime
   pi ~ dbeta(rprime,sprime)
   RplusS <- R+S
   n ~ dbin(pi,RplusS)
   RSn <- (R+S-n)*step(R+S-n)
   RScrs <- logfact(n)+logfact(RSn)-logfact(R+S)
   for (j in 1:N){
       Rj[j] <- (R-j+1)*step(R-j+1)
       Rcr[j] <- logfact(j-1)+logfact(Rj[j])-logfact(R)
       nj[j] <- (n-j+1)*step(n-j+1)
       Sj[j] <- (S-n+j-1)*step(S-n+j-1)
       Scs[j] <- logfact(nj[j])+logfact(Sj[j])-logfact(S)
       pp[j] <- Rcr[j]+Scs[j]-RScrs
       qq[j] <- exp(-pp[j])*step(n+1-j)
    }
   for (j in 1:N){
      phyp[j] <- qq[j]/sum(qq[])
   }
   r ~ dcat(phyp[])
   s <- n-r
}

data;
list(N=500,R=41,n=32,r=8,Rprime=25,rprime=2)