mothMortality.Bayes = function() {

# -------------data------------------------------------
  N = 20
  y = c(1,4,9,13,18,20, 0,2,6,10,12,16)
  sex = c(rep('male',6), rep('female',6))
  dose = rep(c(1,2,4,8,16,32), 2)

  ldose = log(dose)/log(2)
  sexcode = rep(0,length(sex))
  i = sex=='male'
  sexcode[i] = 1

# -------------arguments for R2WinBUGS---------------------
  data = list(n=length(y), N=N, y=y, x=ldose, z=sexcode)
  params = list('alpha','beta','gamma','w')
  inits = function() {
    list(alpha=rnorm(1), beta=rnorm(1), gamma=rnorm(1), w=rbinom(1,1,0.5))
  }

# -------------native WinBUGS code---------------------------
modelFilename = 'model.txt'
cat('
model {

alpha ~ dnorm(0, 0.01)
beta ~ dnorm(0, 0.01)
w ~ dbin(0.5,1)
#tausq <- 0.01
#prec <- w*tausq + (1-w)*10/tausq
#gamma ~ dnorm(0, prec)
gamma ~ dnorm(0, 0.01)

for (i in 1:n) {
    y[i] ~ dbin(p[i], N)
    logit(p[i]) <- alpha + beta*x[i] + w*gamma*z[i]
}
}
', fill=TRUE, file=modelFilename)

# ------------call bugs() to fit model-----------------------
  library(R2WinBUGS)
  fit = bugs(data, inits, params, model.file=modelFilename,
    n.chains=1, n.iter=100000, n.burnin=50000, n.thin=5,  bugs.seed=sample(1:9999,size=1), debug=FALSE, DIC=FALSE)

  post = fit$sims.matrix
  post
}