N <- 100000
epsilon <- 3
sdev <- 0.05
k <- 0
etaset <- 0.5
xstar <- c(125,18+20,34)
n <- sum(xstar)
eta <- rep(NA,N)
for (j in 1:N) {
  etatrial = rnorm(1,etaset,sdev)
  if ((etatrial>0)&(etatrial<1)) {
    inds <- sample(3,n,replace=T,
              p=c((0.5+etatrial/4),((1-etatrial)/2),
              (etatrial/4)))
    x <- c(sum(inds==1),sum(inds==2),sum(inds==3))
    dist <- sqrt(sum((xstar-x)^2))
    if (dist <= epsilon) {
      etaset <- etatrial
      k <- k + 1
    }
  }
eta[j] <- etaset
}
m <- mean(eta)
s <- sd(eta)
cat("k",k,"m",m,"s",s,"\n")