model { log.alpha ~ dnorm( 0.0, 1.0E-6 ) log.beta ~ dnorm( 0.0, 1.0E-6 ) for ( i in 1:n ) { y[ i ] ~ dgamma( alpha, beta ) } alpha <- exp( log.alpha ) beta <- exp( log.beta ) theta <- alpha / beta y.new ~ dgamma( alpha, beta ) }