set.seed(1)
s <- 1e5
K <- seq (5, 100, 5)
original <- augmented <- matrix(0, length(K), 1)
rownames(original) <- rownames(augmented) <- K
row <- 0
for (k in K) {
row <- row + 1
for (sim in 1:s) {
#Generate data from a pupil on k items
theta <- rnorm (1)
b <- runif (k, -1,1)
x <- sum (rlogis (k) <= theta - b)
#Generate a proposal distribution
z = rlogis (k + 1, location = c(b, 0), scale = 1)
j = x + 1
o = order(z)
theta.ast = z[o[j]] #choose (x + 1)-largest as proposed point
#Compute log of acceptance probability
alpha = 0
if(o[j] <= k) {
alpha = alpha + log ( 1 + exp (theta.ast - b[o[j]]))
alpha = alpha + log ( 1 + exp (theta))
alpha = alpha - log ( 1 + exp (theta - b[o[j]]))
alpha = alpha - log ( 1 + exp (theta.ast))
}
if (log (runif (1)) <= alpha)
augmented [row] <- augmented [row] + 1/s
#Original SVE algorithm
theta.ast <- rnorm (1)
x.ast <- sum (rlogis (k) <= theta.ast - b)
if (log (runif (1)) <= (theta - theta.ast)*(x.ast - x))
original [row] <- original [row] + 1/s
}
}
#Figure 6
par (mfrow=c(1,2), cex.main = 1.5, mar = c(4, 5, 0, 0) + 0.1, mgp = c(3.5, 1, 0),
cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot (K, original, ylim = c(0,1) , xlim = c(0,100), lwd = 2, lty = 1,
xlab = " ", ylab = " ", pch = 20)
mtext ("Number of items",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,100,0.2,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.4,100,0.4,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.6,100,0.6,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.8,100,0.8,lty = 2, col = "gray", lwd = 0.5)
segments (0,1,100,1,lty = 2, col = "gray", lwd = 0.5)
plot(K, augmented, ylim = c(0,1) , xlim = c(0,100), lwd = 2, lty = 1, xlab = " ",
ylab = " ", pch = 20)
mtext ("Number of items",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,100,0.2,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.4,100,0.4,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.6,100,0.6,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.8,100,0.8,lty = 2, col = "gray", lwd = 0.5)
segments (0,1,100,1,lty = 2, col = "gray", lwd = 0.5)