#This file contains all the proposed procedures applied to the Rasch # model. It is assumed here that the n pupil abilities have a standard # normal prior, i.e., pi(theta) = normal(0, 1). The k item difficulties # are sampled from a uniform distribution. #Generating a test score for a single pupil on k items k <- 20 #items theta <- rnorm (1) delta <- runif (k, -1, 1) x <- sum (rlogis (k) <= theta - delta) #Generating test scores for n pupils on k items n <- 1e5 #pupils k <- 20 #items theta <- rnorm (n) #pupil abilities delta <- runif (k, -1, 1) #item difficulties x <- rep (0, n) for (i in 1:k) x <- x + 1*(rlogis(n) <= theta - delta[i]) #Original SVE algorithm for a single pupil ## Current state is theta.prime ## Proposal theta.ast <- rnorm (1) #normal prior x.ast <- sum (rlogis (k) <= theta.ast - delta) ## Metropolis step log.p.accept <- (theta.ast - theta.prime) * (x - x.ast) if (log (runif (1)) <= log.p.accept) theta.prime <- theta.ast ## Return theta.prime #Oversampling SVE algorithm for a single pupil ## Current state is theta.prime ## Number of proposals m <- 5 ## Proposal theta.ast <- rnorm (m) #normal prior x.ast <- rep (0, m) for (i in 1:m) x.ast [i] <- sum (rlogis (k) <= theta.ast[i] - delta) i <- which (abs(x.ast - x) == min (abs (x.ast-x))) x.ast <- x.ast [i] theta.ast <- theta.ast [i] ## Metropolis step log.p.accept <- (theta.ast - theta.prime) * (x - x.ast) if (log (runif (1)) <= log.p.accept) theta.prime <- theta.ast ## Return theta.prime #Matching SVE algorithm for n pupils ## Current state is theta.prime ## Proposal theta.ast <- rnorm (n) #normal prior x.ast <- rep (0, n) for (i in 1:k) x.ast <- x.ast + 1*sum (rlogis (n) <= theta.ast - delta[i]) ## Order proposals and targets w.r.t. x and x.ast O <- order(order (x)) #this needs to be done only once! o <- order(x.ast) theta.ast <- theta.ast [O[o]] x.ast <- x.ast[O[o]] ## Metropolis step: log.p.accept <- (theta.ast - theta.prime) * (x - x.ast) log.u <- log(runif (n)) theta.prime [log.p.accept > log.u] <- theta.ast [log.p.accept > log.u] ## Return theta.prime (vector length n) #Rejection algorithm for a single pupil ## Rejection repeat { theta <- rnorm (1) #normal prior ax <- sum (rlogis (k) <= theta - delta) if (ax == x) break } ## Return theta #Rejection algorithm with recycling for n pupils ## Precompute (needs to be done only once) nX <- sapply (0:k, function (score) sum (x == score)) ## Rejection theta.ast <- x.ast <- vector (length = n) while (sum (nX) > 0) { atheta <- rnorm (1) #normal prior ax <- sum (rlogis (k) <= atheta - delta) if (nX [ax + 1] > 0) { nX [ax + 1] <- nX [ax + 1] - 1 theta.ast [sum (nX) + 1] <- atheta x.ast [sum (nX) + 1] <- ax } } ## Assigning generated values (if necessary) for (score in 0:k) theta.prime [x == score] <- theta.ast [x.ast == score] ## Return theta.prime (vector length n) #Binning SVE algorithm for one pupil ## Current state is theta.prime ## Bin size a <- 1 ## Rejection step for proposal repeat { theta.ast <- rnorm (1) #normal prior x.ast <- sum (rlogis (k) <= theta.ast - delta) if (abs (x.ast - x) <= a) break } ## Metropolis step log.p.accept <- (theta.ast - theta.prime) * (x - x.ast) if (log (runif (1)) <= log.p.accept) theta.prime <- theta.ast ## Return theta.prime #Data augmented SVE algorithm for one pupil ## Current state is theta.prime ## Proposal z = rlogis (k + 1, location = c(delta, 0), scale = 1) j = x + 1 o = order(z) theta.ast = z[o[j]] #choose (x + 1)-largest as proposed point ## Metropolis step log.p.accept <- 0 if(o[j] <= k) { log.p.accept = log.p.accept + log (1 + exp (theta.ast - delta[o[j]])) log.p.accept = log.p.accept + log ( 1 + exp (theta.prime)) log.p.accept = log.p.accept - log ( 1 + exp (theta.prime - delta[o[j]])) log.p.accept = log.p.accept - log ( 1 + exp (theta.ast)) } if (log (runif (1)) <= log.p.accept) theta.prime <- theta.ast ## Return theta.prime