N <- 100000
etastar <- c(125,18+20,34)
n <- sum(etastar)
eps <- 3
d <- rep(NA,N)
trial <- rep(NA,N)
for (j in 1:N) {
  etatrial <- runif(1)
  inds <- sample(3,n,replace=T,
     p=c(0.5+etatrial/4,(1-etatrial)/2,etatrial/4))
  samp <- c(sum(inds==1),sum(inds==2),sum(inds==3))
  d <- sqrt(sum((etastar-samp)^2))
  if (d <= eps) trial[j] <- etatrial
}
eta <- trial[!is.na(trial)]
k <- length(eta)
m <- mean(eta)
s <- sd(eta)
cat("k",k,"m",m,"s",s,"\n")