set.seed(1)

s <- 1e4 #number of simulations (see below)
N <- seq(25, 1000, 25) #numbers of pupils used in simulations
k <- 20 #number of items
b <- runif (k, -1, 1) #item difficulties

original <- matching <- matrix (0, length (N), 1)
rownames (original) <- rownames (matching) <- N
row <- 0 #row counter

for (n in N) { #loop over number of pupils N
	accept.sve <- accept.matching <- 0
	for (sim in 1:s) {
		theta <- rnorm (n) #true abilities of n pupils
		theta.ast <- rnorm (n) #proposed states
		x <- x.ast <- rep (0, n)
		for(i in 1:k) {
			#``Observed'' test scores of n pupils
			 x <- x + 1*(rlogis (n) <= (theta - b[i])) 
			#Test scores for the n proposals
			x.ast <- x.ast + 1*(rlogis (n) <= (theta.ast - b[i])) 
		}
		#SVE acceptance rate:
		accept.sve <- accept.sve + 
		  sum (log (runif (n)) <= (x.ast - x)*(theta - theta.ast))

		#Matching proposals to targets (order test scores)
		O <- order(order(x))
		o <- order(x.ast)
		theta.ast <- theta.ast[o[O]]
		x.ast <- x.ast[o[O]]

		#SVE, with matching, acceptance rate
		accept.matching <- accept.matching + 
		  sum (log (runif (n)) <= (x.ast - x)*(theta - theta.ast))
	}
	row <- row + 1
	original [row] <- accept.sve / (n * s)
	matching [row] <- accept.matching / (n * s)
}

#Figure 3
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, matching, ylim = c(0,1) , 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 ("Acceptance rate", 2, line = 3.5, cex = 1.5, font = 2, las = 0)
segments (0,0.2,1000,0.2,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.4,1000,0.4,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.6,1000,0.6,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.8,1000,0.8,lty = 2, col = "gray", lwd = 0.5)
segments (0,1,1000,1,lty = 2, col = "gray", lwd = 0.5)

clip (0,1000,0,1)
abline (h = mean (original), lwd = 4)