ovenCounts.Bayes = function() { # read data from file ymat = as.matrix(read.table('oven.txt')) # model specification in WinBUGS modelFilename = "ovenBayes.txt" cat(" model { beta ~ dnorm(0, 1.0E-2) lambda <- exp(beta) p ~ dunif(0,1) for (i in 1:n) { N[i] ~ dpois(lambda) for (k in 1:K) { y[i,k] ~ dbin(p, N[i]) } } } ", fill=TRUE, file=modelFilename) # arguments for bugs() data = list(n=dim(ymat)[1], K=dim(ymat)[2], y=ymat) params = list('lambda', 'p') inits = function() { ySum = apply(ymat, 1, max) lambdaGuess = mean(ySum) list(p=runif(1,0,1), beta=rnorm(1,mean=log(lambdaGuess),sd=1), N=ySum+1) } # call to bugs() library(R2WinBUGS) fit = bugs(data, inits, params, model.file=modelFilename, n.chains=1, n.iter=30000, n.burnin=10000, n.thin=5, bugs.seed=sample(1:9999,size=1), debug=FALSE, DIC=FALSE) fit }