########################################################################
#This program calculates ability estimates given known item parameters #
#using MAP                                                             #
########################################################################

#for this example, there are just two items, whose parameters are below:
a <- c(1,2)	# a parameters
b <- c(0,1)	# b parameters


thetaList <- seq(-3,3,0.1)	# list of theta values for plots and quadrature

Tlist <- list(0,0)			# lists of values of the trace lines

#Loops through all the items.
for (i in 1:length(a)) {

#compute probabilities of a correct responde using the 2 PL 
#given the theta parameters in thetaList
#and the item parametrs in the a and b vectors.
	Tlist[[i]] <- 1/(1+exp(-a[i]*(thetaList-b[i])))
}

# draw something like Test Scoring's Figure 3.10, in three panels:
#quartz(width=8, height=10)	# note: quartz is specific to Mac OS X; omit on PCs
par(mfrow=c(3,1), cex.axis=2, cex.lab=2, mar=(c(5, 8, 4, 2)+0.1), lab=c(7,3,7))

#Plot the ICC of the first item.
plot(thetaList,Tlist[[1]],type="l", xlab="Theta", ylab="T(u=1)", col="blue")

#Plot the ICC of the second item.
plot(thetaList,(1-Tlist[[2]]),type="l", xlab="Theta", ylab="T(u=0)", col="blue")

#compute the likelihood = P(1-P)
likelihood <- Tlist[[1]]*(1-Tlist[[2]])

#plots the likelihood as a function of theta.
plot(thetaList,likelihood,type="l", xlab="Theta", ylab="Likelihood", col="blue")

################################################################################

#vector of responses to the two items for a person.
#this person responded the first item correctly and the second incorrectly.
#for MLE, you need a minimum of 2 responses, one correct and one incorrect.
u <- c(1,0)	# responses, positive (right) and negative (wrong)

# now draw the log of the likelihood:

loglikelihood <- log(likelihood)
plot(thetaList,loglikelihood,type="l", xlab="Theta", ylab="loglikelihood", col="blue")


################################################################################
 # MAXIMUM A POSTERIORI (MAP) ESTIMATION
 
#Sets graphing mode
#quartz(width=8, height=10)
par(cex.axis=2, cex.lab=2, mar=(c(5, 8, 4, 2)+0.1), lab=c(7,3,7))

#define a population distribution (prior distribution)
popdist <- exp(-0.5*(thetaList^2))	# compute (proportional to) pop. dist.

# re-compute and re-draw the log of the likelihood with the population distribution.
loglikelihood <- log(likelihood) + log(popdist)
plot(thetaList,loglikelihood,type="l", xlab="Theta", ylab="loglikelihood", col="blue")

# Use the "fixed" version of the Newton-Raphson algorithm with stepsize limit and backtracking:
thetahat <- 0	# thetahat starts at 0, all else same as above
Niter <- 10
StepLimit <- 0.5	# this will limit change to 0.5 per iteration
Lastdl <- 0			# and we need to know last direction (start with none)

#loop through iterations for MAP estimation.
for (iter in 1:Niter) {
	l <- 0		# initialize loglikelihood (don't need to compute this,
	dl <- 0		# except for the graphics), and first and second derivatives
	d2l <- 0
	
	
#loop through responses
	for (i in 1:length(u)) {
	
	#obtain the probability of correct response to an item given the 2PL and item parameters.
		T <- 1/(1+exp(-a[i]*(thetahat-b[i])))
		
		#l is the loglikelihood of the response pattern.
	#it is obtained with a summation of l for a response with the l of previous responses.
		l <- l + u[i]*log(T) + (1-u[i])*log(1 - T)
		
	#dl is the first partial derivative of the loglikelihood with respect to theta
  #it is obtained with a summation of dl of a response to the dl of previous reponses 		
		dl <- dl + a[i] * (u[i]-T)  
		
		 #d2l is the second partial derivative of the loglikelihood with respect to theta.
 #it is obtained with a summation (of negative values) of d2l of a response to the d2l of previous responses.
		d2l <- d2l - (a[i]^2) * T * (1 - T)
	}

# add population distribution terms into loglikelihood, first and second derivatives
	l <- l - 0.5*(T^2) # pop. dist.
	dl <- dl - thetahat
	d2l <- d2l - 1

  #calculate the amount to update theta.
	correction <- dl/d2l

#========================  		  
# the twelve lines below are purely for graphics
	if((abs(correction) < 0.000001) | (iter == Niter))
		lcolor <- "red" else lcolor <- "pink"
	ltheta <- thetahat - 0.5	# these lines draw the
	rtheta <- thetahat + 0.5	# straight lines that
	llinear <- l - dl*0.5		# illustrate the gradient
	rlinear <- l + dl*0.5
	lines(c(ltheta,rtheta),c(llinear, rlinear), col=lcolor)
# the five lines below compute and plot the "fitted" quadratic
	q2 <- d2l/2
	q1 <- dl - (2*q2*thetahat)
	q0 <- l - q1*thetahat - q2*thetahat*thetahat
	quad <- q0 + q1*thetaList + q2*thetaList*thetaList
	lines(thetaList, quad, col=lcolor)
 #============================================================
 
	print(c(thetahat, l, dl, d2l))

	if (d2l > -0.01) d2l <- -0.01 # insurance: don't (get close to) zero divide
 
 # below impose stepsize limit
	if (abs(correction) > StepLimit) correction <- sign(correction)*StepLimit
	
	#updates theta.
	thetahat <- thetahat - correction
	
	#Compare the size of the update amount with the convergence criterion.
	if(abs(correction) < 0.000001) break
	
	if ((dl*Lastdl) < 0.0) {					# if gradient has changed sign
		thetahat <- thetahat + (0.5*correction)	# back up half way
		StepLimit <- 0.5*StepLimit				# "split the difference" or
		dl <- 0									# "false position"
	 }
	Lastdl <- dl
}