set.seed(1)
n <- 1e6 #number of pupils
k <- 20 #number of items
#Generate test scores (x) for n pupils on k items
b <- runif (k, -1, 1) #item difficulty
t <- rnorm (n) #normal prior or distribution pi(theta)
x <- rep (0, n) #the test scores X_+
for (i in 1:k)
x <- x + 1*(rlogis (n) <= (t - b[i]))
#Compute the mixture probabilities
px <- sapply (0:k, function(s) sum (x == s))/n
#Figure 1
par(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(0:k, px,
ylim = c(0,1) , xlim = c(0,20), lwd = 2, lty = 1, xlab = " ",
ylab = " ", pch = 20)
mtext(expression(x["+"]),1, line = 3, cex = 1.5, font = 2)
mtext(expression(paste(pi,"(",x["+"],")")), 2, line = 3.5, cex = 1.5,
font = 2, las = 0)
segments (0,0.2,20,0.2,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.4,20,0.4,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.6,20,0.6,lty = 2, col = "gray", lwd = 0.5)
segments (0,0.8,20,0.8,lty = 2, col = "gray", lwd = 0.5)
#compute acceptance rate for pupil with test score 9
theta <- 0
rate <- 0
for (i in 1:1e5) {
theta.ast <- rnorm (1)
x.ast <- sum (rlogis(k) <= theta.ast - b)
if (log (runif (1)) <= (theta - theta.ast) * (x.ast - 9)) {
theta <- theta.ast
rate <- rate + 1
}
}
rate / i
#Probability to generate kernel with |tx.ast - 9| > 4
1 - sum (px[6:14])