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)