#######################################################################
#This program calculates ability estimates given known item parameters#
#using MLE.                                                           #
#######################################################################

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

# create a vector of equally spaced theta values
thetaList <- seq(-3,3,0.1)

# create an empty list to store probabilities.
Problist <- list(0,0)

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

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

} #close loop through items

#set graphic options to display three graphs at once.
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,Problist[[1]],type="l", xlab="Theta", ylab="T(u=1)", col="blue")

#Plot the ICC of the second item.
plot(thetaList,Problist[[2]],type="l", xlab="Theta", ylab="T(u=0)", 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)

#compute the likelihood function
likelihood <- (Problist[[1]]^u[1]*(1-Problist[[1]])^(1-u[1]))*
              (Problist[[2]]^u[2]*(1-Problist[[2]])^(1-u[2]))

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

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

# now draw the log of the likelihood:
par(cex.axis=2, cex.lab=2, mar=(c(5, 8, 4, 2)+0.1), lab=c(7,3,7)) #set graphic parameters
loglikelihood <- log(likelihood)
plot(thetaList,loglikelihood,type="l", xlab="Theta", ylab="loglikelihood", col="blue")

#========================
#MAXIMUM LIKELIHOOD ESTIMATION for a person given the responses for the two items defined
#previously 
# set up for  Newton-Raphson iterations:
thetahat <- 0	# thetahat starts at zero
# iterations begin below, maximum of 10 often enough
Niter <- 10

#loops through iterations for Maximum likelihood 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
	
	#loops through items (which is the same as looping through responses).
	for (i in 1:length(a)) {
	
	#T is the probability function for an item with the 1PL or 2PL
		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)
	
  #finishes loop through responses for a person.
  }

  #---------------------
# the twelve lines below are purely for graphics
	if((abs(dl) < 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)
#----------------------------

#provides the current estimate of theta, the loglikelihood, the first and secon partial derivatives.
	print(c(thetahat, l, dl, d2l))
	
#terminates the loop (MLE converges) if  the first derivative is smaller than the criterion.
	if(abs(dl) < 0.000001) break	

  #if the convergence criterion above is not met, updates the estimate of thetahat
  #with the ration of first and second partial derivatives.
	thetahat <- thetahat - (dl/d2l)

#FINISHES MLE ESTIMATION.
}