#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