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

#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 the loglikelihood
par(mfrow=c(1,1), cex.axis=2, cex.lab=2, mar=(c(5, 8, 4, 2)+0.1), lab=c(7,3,7))
plot(thetaList,loglikelihood,type="l", xlab="Theta", ylab="loglikelihood", col="blue")


################################################################################
#EAP ESTIMATION

# graphing parameters.
#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))

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

#normalize the density of the population distribution (making it sum to 1)
popdist <- popdist / sum(popdist)

#plot population distribution.
plot(thetaList,popdist,type="l", xlab="Theta", ylab="Pop. Dist.", col="blue")

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

#add ICC of item 2.
lines(thetaList,(1-Tlist[[2]]),type="l", col="red")

#compute likelihood adjusted by the density of the population distribution. 
likelihood <- Tlist[[1]]*(1-Tlist[[2]])*popdist

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

#------------------------
# add quadrature histobars to graphic:

# compute difference between quadrature points
dtheta <- thetaList[2] - thetaList[1]	

halfdtheta <- dtheta/2

#loop through abilities values.
for (i in 1:length(thetaList)) {

#add rectangles to the histogram for each quadrature point.
	lines(c(thetaList[i]-halfdtheta, thetaList[i]+halfdtheta), 
			c(likelihood[i], likelihood[i]), col="magenta")
	lines(c(thetaList[i]-halfdtheta, thetaList[i]-halfdtheta),
			c(likelihood[i], 0), col="magenta")
	lines(c(thetaList[i]+halfdtheta, thetaList[i]+halfdtheta),
			c(likelihood[i], 0), col="magenta")
}

#-------------------------------------
# compute EAP[theta] and SD[theta]
EAPtheta <- sum(likelihood*thetaList) / sum(likelihood)


SDtheta <- sqrt(sum(likelihood*((thetaList - EAPtheta)^2)) / sum(likelihood))
print(c(EAPtheta, SDtheta))
# add to the graphic
MaxLikelihood <- max(likelihood)
lines(c(EAPtheta, EAPtheta), c(MaxLikelihood, 0), col="black")
text(EAPtheta, 0, round(100*EAPtheta)/100, pos=4, cex=1.5, col="black")
lines(c(EAPtheta, EAPtheta+SDtheta), c(0.6*MaxLikelihood, 0.6*MaxLikelihood),
				 col="black")
lines(c(EAPtheta, EAPtheta-SDtheta), c(0.6*MaxLikelihood, 0.6*MaxLikelihood),
				 col="black")
text(EAPtheta+0.5*SDtheta, 0.6*MaxLikelihood, round(100*SDtheta)/100, 
			pos=3, cex=1.5, col="black")

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

# re-do it with one-fifth that many quadrature points
thetaList <- seq(-3,3,0.5)	# list of theta values for plots and quadrature

Tlist <- list(0,0)			# lists of values of the trace lines
for (i in 1:length(a)) {
	Tlist[[i]] <- 1/(1+exp(-a[i]*(thetaList-b[i])))
}
popdist <- exp(-0.5*(thetaList^2))	# compute (proportional to) pop. dist.

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

popdist <- popdist / sum(popdist) # normalize pop. dist, for graphics and later
plot(thetaList,popdist,type="l", xlab="Theta", ylab="Pop. Dist.", col="blue")

plot(thetaList,Tlist[[1]],type="l", xlab="Theta", ylab="T(u=1)", col="red")
lines(thetaList,(1-Tlist[[2]]),type="l", col="red") 
likelihood <- Tlist[[1]]*(1-Tlist[[2]])*popdist
plot(thetaList,likelihood,type="l", xlab="Theta", ylab="Likelihood", col="magenta")

# add quadrature histobars to graphic:
# compute difference between quadrature points
dtheta <- thetaList[2] - thetaList[1]	
halfdtheta <- dtheta/2
for (i in 1:length(thetaList)) {
	lines(c(thetaList[i]-halfdtheta, thetaList[i]+halfdtheta), 
			c(likelihood[i], likelihood[i]), col="magenta")
	lines(c(thetaList[i]-halfdtheta, thetaList[i]-halfdtheta),
			c(likelihood[i], 0), col="magenta")
	lines(c(thetaList[i]+halfdtheta, thetaList[i]+halfdtheta),
			c(likelihood[i], 0), col="magenta")
}

# compute EAP[theta] and SD[theta]
EAPtheta <- sum(likelihood*thetaList) / sum(likelihood)
SDtheta <- sqrt(sum(likelihood*((thetaList - EAPtheta)^2)) / sum(likelihood))
print(c(EAPtheta, SDtheta))


# add to the graphic
MaxLikelihood <- max(likelihood)
lines(c(EAPtheta, EAPtheta), c(MaxLikelihood, 0), col="black")
text(EAPtheta, 0, round(100*EAPtheta)/100, pos=4, cex=1.5, col="black")
lines(c(EAPtheta, EAPtheta+SDtheta), c(0.6*MaxLikelihood, 0.6*MaxLikelihood),
				 col="black")
lines(c(EAPtheta, EAPtheta-SDtheta), c(0.6*MaxLikelihood, 0.6*MaxLikelihood),
				 col="black")
text(EAPtheta+0.5*SDtheta, 0.6*MaxLikelihood, round(100*SDtheta)/100, 
			pos=3, cex=1.5, col="black")