set.seed (1)

s <- 1e4
N <- seq (25, 1000, 25)
k <- 20
b <- runif (k, -1,1) #item difficulties

#The original rejection algorithm's efficiency does not depend on n, 
# so we will do a single long run to estimate its efficiency.

n <- max (N)
count <- 0
for (sim in 1:s) {
	#Generate item responses from n pupils on k items
	theta <- rnorm (n)
	x <- rep (0, n)
	for(i in 1:k) 
		x <- x + 1*(rlogis(n) <= (theta - b[i]))
	
	for (p in 1:n) {
		x.ast <- -1
		while (x.ast != x[p]) {
			count <- count + 1
			theta.ast <- rnorm (1)
			x.ast <- sum (rlogis (k) <= (theta.ast - b))
		}
	}
}
rejection <- count / (n * s)

#Recycling:
recycling <- matrix (0, length (N), 1)
row <- 0
for (n in N) {
	count <- 0
	row <- row + 1
	for (sim in 1:s) {
		#Generate item responses from n pupils on k items
		theta <- rnorm (n)
		x <- rep (0, n)
		for(i in 1:k) 
			x <- x + 1*(rlogis(n) <= (theta - b[i]))
		nX <- sapply (0:k, function (score) sum (x == score))
		
		while (sum(nX) > 0) {
			count <- count + 1
			theta.ast <- rnorm (1)
			x.ast <- sum (rlogis (k) <= (theta.ast - b))
			if (nX [x.ast + 1] > 0)
				nX [x.ast + 1] <- nX [x.ast + 1] - 1
		}
	print (paste (sim/s, "    ", row / length (N)))
	}
	recycling[row] <- count / (n * s)
}

#Figure 4
par (cex.main = 1.5, mar = c(4, 5, 0, 0.5) + 0.1, mgp = c(3.5, 1, 0), 
  cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(N, recycling - 1, ylim = c(0,25) , xlim = c(0,1000), lwd = 2, lty = 1, 
  xlab = " ", ylab = " ", pch = 20)
mtext("Number of latent variables",1, line = 3, cex = 1.5, font = 2)
mtext("Number of rejected samples", 2, line = 3.5, cex = 1.5, font = 2, las = 0)
for (i in 1:5) 
	segments (0,i*5,1000,i*5,lty = 2, col = "gray", lwd = 0.25)

par (new = TRUE)
points (N, recycling - 1, pch = 19)
clip (0,1000,0,25)
abline (h = rejection, lwd = 4)