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)