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
for (i in 1:length(a)) {
	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(thetaList,Tlist[[1]],type="l", xlab="Theta", ylab="T(u=1)", col="blue")
plot(thetaList,(1-Tlist[[2]]),type="l", xlab="Theta", ylab="T(u=0)", col="blue")
likelihood <- Tlist[[1]]*(1-Tlist[[2]])
plot(thetaList,likelihood,type="l", xlab="Theta", ylab="Likelihood", col="blue")

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

u <- c(1,0)	# responses, positive (right) and negative (wrong)

# now draw the log of the likelihood:
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))
loglikelihood <- log(likelihood)
plot(thetaList,loglikelihood,type="l", xlab="Theta", ylab="loglikelihood", col="blue")


# set up for (un"fixed") Newton-Raphson iterations:
thetahat <- 0	# thetahat starts at zero
# iterations begin below, maximum of 10 often enough
Niter <- 10
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
	for (i in 1:length(a)) {
		T <- 1/(1+exp(-a[i]*(thetahat-b[i])))
		l <- l + u[i]*log(T) + (1-u[i])*log(1 - T)
		dl <- dl + a[i] * (u[i]-T)  
		d2l <- d2l - (a[i]^2) * T * (1 - T)
	}

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

	print(c(thetahat, l, dl, d2l))

	if(abs(dl) < 0.000001) break	#call
	thetahat <- thetahat - (dl/d2l)
}


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

# now re-draw the log of the likelihood, with the plot re-scaled 
# to leave room for grief:
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))
loglikelihood <- log(likelihood)
plot(thetaList,loglikelihood,type="l", xlab="Theta", ylab="loglikelihood", 
		ylim=c(-4,2), col="blue")

# now show it "breaking" by starting at -1.25:
thetahat <- -1.25	# thetahat starts at -1.25, all else same as above
Niter <- 5
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
	for (i in 1:length(a)) {
		T <- 1/(1+exp(-a[i]*(thetahat-b[i])))
		l <- l + u[i]*log(T) + (1-u[i])*log(1 - T)
		dl <- dl + a[i] * (u[i]-T)  
		d2l <- d2l - (a[i]^2) * T * (1 - T)
	}

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

	print(c(thetahat, l, dl, d2l))

	if(abs(dl) < 0.000001) break	#call
	thetahat <- thetahat - (dl/d2l)
}

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

# re-draw the log of the likelihood: 
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))
loglikelihood <- log(likelihood)
plot(thetaList,loglikelihood,type="l", xlab="Theta", ylab="loglikelihood", 
		ylim=c(-4,2), col="blue")

# now "fix" it, with stepsize limit and backtracking starting at -1.25:
thetahat <- -1.25	# thetahat starts at -1.25, 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)
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
	for (i in 1:length(a)) {
		T <- 1/(1+exp(-a[i]*(thetahat-b[i])))
		l <- l + u[i]*log(T) + (1-u[i])*log(1 - T)
		dl <- dl + a[i] * (u[i]-T)  
		d2l <- d2l - (a[i]^2) * T * (1 - T)
	}

	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
	thetahat <- thetahat - correction
	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
}

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

# re-compute and re-draw the log of the likelihood with the population dist.:
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))
popdist <- exp(-0.5*(thetaList^2))	# compute (proportional to) pop. dist.
loglikelihood <- log(likelihood) + log(popdist)
plot(thetaList,loglikelihood,type="l", xlab="Theta", ylab="loglikelihood", col="blue")

# again the "fixed" version with stepsize limit and backtracking:
thetahat <- 0	# thetahat starts at -1.25, 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)
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
	for (i in 1:length(a)) {
		T <- 1/(1+exp(-a[i]*(thetahat-b[i])))
		l <- l + u[i]*log(T) + (1-u[i])*log(1 - T)
		dl <- dl + a[i] * (u[i]-T)  
		d2l <- d2l - (a[i]^2) * T * (1 - T)
	}
# add pop dist. terms into first and second derivatives
	l <- l - 0.5*(T^2) # pop. dist.
	dl <- dl - thetahat
	d2l <- d2l - 1

	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
	thetahat <- thetahat - correction
	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
}

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

# now draw something like Test Scoring's Figure 3.15, 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))

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

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

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