# Albert IRT MCMC 
# using the Stouffer-Toby (1951) data 

# The data table:
x1 <- c(rep(T,4),F,rep(T,3),rep(F,3),T,rep(F,4))
x2 <- c(rep(T,3),F,T,T,F,F,T,T,F,F,T,rep(F,3)) 
x3 <- c(T,T,F,T,T,F,T,F,T,F,T,F,F,T,F,F)
x4 <- c(T,F,T,T,T,F,F,T,F,T,T,F,F,F,T,F)
x <- cbind(x1,x2,x3,x4)
r <- c(42,23,6,6,1,24,25,7,4,2,1,38,9,6,2,20)

nitems <- 4    	
n <- sum(r)    	

stouffer <- data.frame(I(x),r)

# ML estimates from BL-StoufferTobyNO-NumH.R are
#p <- c(0.64741, -0.96905, 0.81028, -0.05289,
#	0.94043, -0.00906, 1.21505, 0.76743)
#a <- rep(0, ncol(x))
#d <- rep(0, ncol(x))
#for (i in 1:ncol(x)) {
#  	a[i] <- p[2*(i-1) + 1]
#  	d[i] <- p[2*i]
#}
a <- rep(1, ncol(x))
d <- rep(0, ncol(x))


Z <- matrix(0, ncol = nitems, nrow = n) 
theta <- matrix(0, ncol = 1, nrow = n)

# start theta with standardized sum score
person <- 0
for(p in 1:length(stouffer$r)) {
	for(freq in 1:stouffer$r[p]) {
		person <- person + 1
		for(i in 1:nitems) {
			if (stouffer$x[p,i] == T) {
				theta[person] <- theta[person] + 1
			}
		}
	}
}
theta <- (theta - mean(theta))/sqrt(var(theta)[1,1])

DrawZ	<-	function(a,d,theta,dataframe) {
# T -> Z ~ N(a*theta-d,1) truncated at the left by 0
# F -> Z ~ N(a*theta-d,1) truncated at the right by 0

randunif <-	matrix(runif(nitems*n), ncol = nitems, nrow = n)

Z1 <- matrix(0, ncol = nitems, nrow = n) 
person <- 0
for(p in 1:length(dataframe$r)) {
	for(freq in 1:dataframe$r[p]) {
		person <- person + 1
		for(i in 1:nitems) {
			mu <- a[i]*theta[person] - d[i]
			BB <- pnorm(-mu)
			if (dataframe$x[p,i] == F) {
				tt <- BB*randunif[person,i]
			}
			else {  
				tt <- (1-BB)*randunif[person,i] + BB
			}
			Z1[person,i] <- qnorm(tt) + mu				
		}
	}
}
	return(Z1)
}

DrawTheta <- function(a,Z) {
	theta1 <- matrix(0, ncol = 1, nrow = n)
	v <- 1 / sum(a^2)
	pvar <-	1/v + 1

	for(i in 1:length(Z[,1])) {
		thetahat <- sum(a*(Z[i,]+d)) 
		mu <-	thetahat / pvar
		theta1[i] <- rnorm(1,mu,sqrt(1/pvar)) 
	}
  return(theta1)
}

DrawItemParams <- function(theta, Z) {
	p <- matrix(0,nrow=2*nitems,ncol=1)
	ones <- matrix(-1,nrow=n,ncol=1)
	X <- matrix(c(ones, theta), nrow=n, ncol=2)
	XtXinv <- solve(t(X)%*%X)
	amat <- chol(XtXinv)
	for(i in 1:nitems) {
		meanIparams <- XtXinv %*% t(X) %*% Z[,i]		
		repeat {
			newIparams <- t(amat) %*% matrix(rnorm(2),nrow=2,ncol=1)+meanIparams
			p[2*(i-1) + 1] <- newIparams[2]
			p[2*i] <- newIparams[1]
			if (newIparams[2] < 0) print(i)
			if (newIparams[2] > 0) break
		}
	}
	return(p)
}


niter <- 5000
hardskip <- 1
nrowResults <- niter
nparams <- 2*length(a)

Results <- matrix (0, nrow=nrowResults, ncol=nparams)


system.time( for(iter in 1:niter) {
	for(skipped in 1:hardskip) {
		Z <- DrawZ(a,d,theta,stouffer)
		theta <- DrawTheta(a,Z)
		p <- DrawItemParams(theta,Z)
		for(i in 1:nitems) {
			a[i] <- p[2*(i-1) + 1]
			d[i] <- p[2*i]
		}
	}
	Results[iter,] <- p
}	)	


colnames(Results) <- c("a1", "d1", "a2", "d2", "a3", "d3", "a4", "d4")

library(coda)
thin <- 1
ResultsMCMC <- mcmc(Results, 1, niter, thin)
plot(ResultsMCMC)

summary(ResultsMCMC)  # fails

autocorr(ResultsMCMC)

# below makes it work
thin <- 5
ResultsMCMC <- mcmc(Results, 1, niter, thin)
summary(ResultsMCMC)   # fails
thin <- 10
ResultsMCMC <- mcmc(Results, 1, niter, thin)
summary(ResultsMCMC)  # fails
thin <- 20
ResultsMCMC <- mcmc(Results, 1, niter, thin)
summary(ResultsMCMC)

thin <- 20
ResultsMCMC <- mcmc(Results, 1000, niter, thin)
summary(ResultsMCMC)



#########################rolled-our-own "trace" plots###################
par(mfrow=c(2,nitems))
for (i in 1:nitems) {
	plot(0:1, 0:1, xlim=c(0,length(Results[,i])), ylim=c(0,3), xlab="sequence",
			ylab="a", type = "n")
	points(seq(1:length(Results[,2*(i-1) + 1])),Results[,2*(i-1) + 1])
}
for (i in 1:nitems) {
	plot(0:1, 0:1, xlim=c(0,length(Results[,i])), ylim=c(-3,3), xlab="sequence",
			ylab="d", type = "n")
	points(seq(1:length(Results[,2*i])),Results[,2*i])
}

#########################rolled-our-own summary###################
amean <- matrix(0,nrow=1,ncol=nitems)
astddev <- matrix(0,nrow=1,ncol=nitems)
dmean <- matrix(0,nrow=1,ncol=nitems)
dstddev <- matrix(0,nrow=1,ncol=nitems)

for (i in 1:nitems) {
	amean[i] <- mean(Results[,2*(i-1) + 1])
	astddev[i] <- sqrt(var(Results[,2*(i-1) + 1]))
	dmean[i] <- mean(Results[,2*i])
	dstddev[i] <- sqrt(var(Results[,2*i]))
}

#########################rolled-our-own histograms###################
par(mfrow=c(2,nitems))
for (i in 1:nitems) {
	hist(Results[,2*(i-1) + 1], main="", xlab=paste("a",i))
}
for (i in 1:nitems) {
	hist(Results[,2*i], main="", xlab=paste("d",i))
}