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)