v <- 10000;

drug<-c("a","a","a","a","a","a","a","a","a","a","d","d","d","d","d","d","d","d","d","d","placebo","placebo","placebo","placebo","placebo","placebo","placebo","placebo","placebo","placebo")
drug<-as.factor(drug)

lbs <- c(6,0,2,8,11,4,13,1,8,0,0,2,3,1,18,4,14,9,1,9,13,10,18,5,23,12,5,16,1,20)

boxplot(lbs~drug)

p <- length(levels(drug))

yBar <- rep(NULL,p)
n    <- rep(NULL,p)
mu   <- rep(0,p*v); dim(mu) <- c(v,p);

S2resid <- 0;
for (i in 1:p) {
  yBar[i] <- mean(lbs[drug==levels(drug)[i]]);
  n[i]    <- sum(drug==levels(drug)[i])
  S2resid <- S2resid + 
	sum((lbs[drug==levels(drug)[i]]-yBar[i])^2)
}
nTotal <- sum(n)

tau <- rgamma(v,(nTotal-p)/2,S2resid/2);
for (i in 1:v) {
  for (j in 1:p) {
    mu[i,j] <- rnorm(1,yBar[j],sqrt(1/(tau[i] * n[j])))
  }
}