N <- 100000
etastar <- c(461,130,161,515)
n <- sum(etastar)
eps <- 22
d <- rep(NA,N)
trial <- rep(NA,N)
for (j in 1:N) {
  etatrial <- runif(1)
  inds <- sample(4,n,replace=T,p=c(0.25+etatrial/4,etatrial/4,
          (1-etatrial)/4,(1-etatrial)/4+0.25))
  samp <- c(sum(inds==1),sum(inds==2),sum(inds==3),sum(inds==4))
  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")