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)