require(rjags)
require(R2jags)
require(devtools)
require(jagstools)

N = 1000
x = rnorm(N,5,1)
mu = exp(1 + 0.5*x)
y = rpois(N, mu)

model = function() {
  for(i in 1:N) {
    y[i] ~ dpois(mu[i])
    log(mu[i]) <- alpha + beta*x[i]
  }
  alpha  ~ dnorm(0,0.001)
  beta  ~ dnorm(0,0.001)
}

ds = list(x=x, y=y, N=N)
params = c('alpha','beta','mu')
nchains = 2
inits = list(data.frame(alpha=-2, beta=0.5), data.frame(alpha=-2, beta=0.5))

jags_model = jags(data=ds, param=params, n.chains=nchains, inits=inits,
                  n.iter=2000, n.burnin=500, model.file=model)
r1 = jagsresults(x=jags_model, params=c('alpha'))
jags_model = jags(data=ds, param=params, n.chains=nchains, inits=inits,
                  n.iter=2000, n.burnin=500, model.file=model)
r2 = jagsresults(x=jags_model, params=c('alpha'))

r1 ; r2 # Small differences (mcmc sampling)


# setting the mcmc seeds
set.seed(100)
jags_model = jags(data=ds, param=params, n.chains=nchains, inits=inits,
                  n.iter=2000, n.burnin=500, model.file=model)
r1 = jagsresults(x=jags_model, params=c('alpha'))

set.seed(100)
jags_model = jags(data=ds, param=params, n.chains=nchains, inits=inits,
                  n.iter=2000, n.burnin=500, model.file=model)
r2 = jagsresults(x=jags_model, params=c('alpha'))

r1 ; r2 #results are now identical