N <- 1000
E <- N/2
x0 <- c(125,18+20,34)
n <- sum(x0)
sdev <- 0.01
pi <- function(x) dunif(x,0,1)
q <- function(eta,eta0) dnorm(eta,eta0,sdev)
# PRC1
epsilon <- c(9,3,1)
T <- length(epsilon)
eta <- matrix(NA,N,T)
etaold <- matrix(NA,N,T)
W <- matrix(NA,N,T)
X <- rep(NA,N)
W[1:N,1] <- 1/N
for (t in 1:T) {
  # PRC2.1
  for (i in 1:N) {
    dist <- max(epsilon) + 1
    while (dist > epsilon[t]) {
      if (t==1) etastarstar <- runif(1,0,1)
        if (t>1) {
          etastar <- 
            sample(eta[,t-1],1,replace=TRUE,prob=W[,t-1])
          etastarstar <- -1
          while (!((0<etastarstar)&(etastarstar<1))) {
            etastarstar <- rnorm(1,etastar,sdev)
          }
        }
        inds <- sample(3,n,replace=T,
                  p=c((0.5+etastarstar/4),((1-etastarstar)/2),
                  (etastarstar/4)))
        xstarstar <- c(sum(inds==1),sum(inds==2),sum(inds==3))
        dist <- sqrt(sum((xstarstar-x0)^2))
    }
    # PRC2.2
    eta[i,t] <- etastarstar
  }
  for (i in 1:N) {
    if (t==1) W[1:N,1] <- 1/N
    if (t>1) {
      S <- sum(W[1:N,1]*q(eta[i,t],eta[1:N,t-1]))
      W[i,t] <- pi(eta[i,t])/S
    }
  }
  ESS <- 1/sum(W[1:N,t]^2)
  if (ESS < E) {
    # PRC3
    etaold <- eta
    for (i in 1:N) eta[i,t] <- 
        sample(eta[,t],1,replace=TRUE,prob=W[,t-1])
    W[1:N,t] <- 1/N
  }
}
plot(density(eta[,1]),lty=3,xlim=c(0.45,0.75),ylim=c(0,15),
  xlab="",ylab="",main="")
par(new=T)
plot(density(eta[,2]),lty=2,xlim=c(0.45,0.75),ylim=c(0,15),
  xlab="",ylab="",main="")
par(new=T)
plot(density(eta[,3]),lty=1,xlim=c(0.45,0.75),ylim=c(0,15),
  xlab="",ylab="",main="")
legend(0.5,12.5,c(expression(paste(epsilon," = 9")),
                  expression(paste(epsilon," = 3")),
                  expression(paste(epsilon," = 1"))),lty=3:1)