# Bock and Lieberman IRT MML 
# using the Stouffer-Toby (1951)

# 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)

theta <- seq(-4,4,0.5)

Gaussian.pts <-	function(mu,sigma,theta) {
 	curve <- exp(-0.5*((theta - mu)/sigma)^2)
 	curve <- curve/sum(curve)
 }

# using numerical derivatives, the only difference between this and normal ogive
# is the trace line points function definition
trace.line.pts <-	function(a,b,theta)	{
 	traceline <- 1/(1+exp(-a*(theta-b)))
 }

ll.2pl.ip <- function(p,testdata,theta) {
	for (i in 1:ncol(testdata$x)) {
  		a[i] <- p[2*(i-1) + 1]
  		b[i] <- p[2*i]
	}

  	itemtrace <- matrix(0,nrow=ncol(testdata$x),ncol=length(theta))

  	for (i in 1:length(a)) {
  		itemtrace[i,] <- trace.line.pts(a[i],b[i],theta)
  	}
  	expected <- rep(0,length(testdata$r))
  	for (i in 1:length(testdata$r)) {
  		posterior <- Gaussian.pts(0,1,theta) 
  		for (item in 1:ncol(testdata$x)) {
  			x <- I(testdata$x[i,item])   
  			if (x)
    			posterior <- posterior*itemtrace[item,]  
  			else
    			posterior <- posterior*(1-itemtrace[item,])
 		}
  		expected[[i]] <- sum(posterior)       
	}
 	l <- (-1)*sum(testdata$r*(log(expected))) 
 }

# distant starting values
a <- c(1,1,1,1)
b <- c(0,0,0,0)
p <- rep(0,2*length(a))
for (i in 1:length(a)) {
  	p[2*(i-1) + 1] <- a[i]
  	p[2*i] <- b[i]
}
system.time(result <- nlm(f=ll.2pl.ip,p=p,hessian=TRUE,
	testdata=stouffer,theta=theta))
result

SE.2PL <- sqrt(diag(solve(result$hessian)))
print(SE.2PL)