# Bock and Lieberman IRT MML 
# 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)

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

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

trace.line.pts <-	function(a,c,theta)	{
 	traceline <- pnorm(a*theta - c)
 }

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

  	itemtrace <- matrix(0,nrow=ncol(testdata$x),ncol=length(theta))
# dPda and dPdc will be the derivatives of the cell Ps w.r.t. a and c,
# for each response-pattern cell
  	dPda <- matrix(0,nrow=length(testdata$r),ncol=length(a))
  	dPdc <- matrix(0,nrow=length(testdata$r),ncol=length(a))
# gr will ultimately be the gradient
  	gr <- rep(0,length(p))

  	for (i in 1:length(a)) {
  		itemtrace[i,] <- trace.line.pts(a[i],c[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)  
# what follows computes Bock & Lieberman's equation 9  		
  		for (item in 1:ncol(testdata$x)) {
  			h <- (1/sqrt(2*pi))*exp((-0.5)*((a[item]*theta - c[item])^2))
  			x <- I(testdata$x[i,item])   
  			if (x) {
    			dPda[i,item] <- sum(theta*h*posterior/itemtrace[item,])
    			dPdc[i,item] <- (-1)*sum(h*posterior/itemtrace[item,])
    		}
  			else {
   			   dPda[i,item] <- (-1)*sum(theta*h*posterior/(1-itemtrace[item,]))
    		   dPdc[i,item] <- sum(h*posterior/(1-itemtrace[item,]))
    		}
  		}

	}
 	l <- (-1)*sum(testdata$r*(log(expected))) 
# the next loop accumulates the gradient as per Bock & Lieberman
# near the bottom of page 183 	
 	for (i in 1:length(a)) {
 		gr[2*(i-1) + 1] <- (-1)*sum((testdata$r/expected)*dPda[,i])
 		gr[2*i] <- (-1)*sum((testdata$r/expected)*dPdc[,i])
 	}
# dPdac assembles the columns of dPda and dPdc into a single matrix
# columns for a1, c1, a2, c2, a3, ...
	dPdac <- matrix(0,nrow=length(testdata$r),ncol=length(p))
 	for (i in 1:length(a)) {
 		dPdac[,2*(i-1) + 1] <- dPda[,i]
 		dPdac[,2*i] <- dPdc[,i]
 	}
# h is a matrix-expression to compute the sums of squares and cross-products
# Bock & Lieberman equation 11, except here expressed in matix operations
 	h <- sum(testdata$r) * (t(dPdac) 
 				%*% diag(1/expected, nrow=length(testdata$r), 
				ncol=length(testdata$r)) %*% dPdac)
	attr(l, "gradient") <- gr
	attr(l, "hessian") <- h  
  return(l)
 }

# above works with my own NR, set to minimize (negative gradient, positive H)
# the gradient alone is fine with nlm, but hessian is not although it's very
# close to the numerical hessian

# Starting values:
a <- c(1,1,1,1)
c <- 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] <- c[i]
}

for (iter in 1:25) {
	update <- ll.NormalOgive.ip(p,stouffer,theta)
# compute the gradient and the hessian (above) and extract them (below)
	g <- attr(update,"gradient")
	H <- attr(update,"hessian")
# compute the next value of the parameters (p) using a Newton-Raphson step
	p <- p - g %*% solve(H)
	print(p)
# call it converged if the norm of the gradient is < 10^-6
	if (sum(g^2) < 0.000001) break
}
SE.NR <- sqrt(diag(solve(H)))
print(SE.NR)