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)