# Bock & Jones Constant Method Saltiness Data:
x <- c(.3, .2, .1, -.1, -.2, -.3)
N <- c(49, 47, 50, 48, 48, 49)
r <- c(48, 38, 31, 13, 3, 2)

# this next bit plots the data, so we know that with which we work
quartz(width=8, height=8)	# note: quartz is specific to Mac OS X; omit on PCs
par(cex.axis=2, cex.lab=2, mar=(c(5, 8, 4, 2)+0.1), lab=c(7,3,7))
plot(x, (r/N), xlab="Difference in Salt Concentration", ylab="P('more salty')", col="blue")

# for R procedures that work from individual data (like glm),
# we need "responses" and "SaltConc" values for each person
# this next section manufactures those from the tabulated data entered above
responses <- c(rep(1,r[1]), rep(0,(N[1] - r[1])),
				rep(1,r[2]), rep(0,(N[2] - r[2])),
				rep(1,r[3]), rep(0,(N[3] - r[3])),
				rep(1,r[4]), rep(0,(N[4] - r[4])),
				rep(1,r[5]), rep(0,(N[5] - r[5])),
				rep(1,r[6]), rep(0,(N[6] - r[6])))
SaltConc <- c(rep(x[1], N[1]), rep(x[2], N[2]), rep(x[3], N[3]), 
				rep(x[4], N[4]), rep(x[5], N[5]), rep(x[6], N[6]))

# R's generalized linear models procedure using iteratively reweighted LS
# (do not confuse this glm with GLM elsewhere that stands for
# General Linear Model---in R, that's done in lm (for Linear Model)
glm.out <- glm(responses~SaltConc, family = binomial(link = "probit"))

# add the fitted line to the graphic
alpha <- glm.out$coefficients[1]
beta <- glm.out$coefficients[2]
lines(x, pnorm(alpha + (beta*x)), col="blue")



# for derivative-free ML (R minimizes the negative loglikelihood):
ll.probit <- function(params, x, N, r) {
	alpha <- params[1]
	beta <- params[2]
	P <- pnorm(alpha + (beta*x))
	l <- (-1)*sum(r*log(P) + ((N-r)*log(1-P)))
 }

params <- c(0,5)

system.time( result <- nlm(f=ll.probit,p=params,hessian=TRUE,x=x,N=N,r=r) )
print(result)
SE <- sqrt(diag(solve(result$hessian)))
print(SE)


# while we're here, note that the ll.probit routine above 
# was used to make the likelihood contour graphics as well, as follows:
### plot likelihood (as contours) for alpha-beta ###
x.alpha <- seq((result$estimate[1]-3.5*SE[1]),
		 (result$estimate[1]+3.5*SE[1]), length=100)
y.beta <- seq((result$estimate[2]-3.5*SE[2]),
		 (result$estimate[2]+3.5*SE[2]), length=100)
z <- matrix(0, nrow=100, ncol=100)
for (i in 1:100) {
	for (j in 1:100) {
		params <- c(x.alpha[i], y.beta[j])
		z[i,j] <- ll.probit(params, x, N, r)
	}
 }

NMLL<- result$minimum
quartz(width=8, height=8)	# note: quartz is specific to Mac OS X; omit on PCs
par(cex.axis=2, cex.lab=2, cex.main=2, mar=(c(5, 8, 4, 2)+0.1), lab=c(7,7,7))
contour(x.alpha, y.beta, z, 
   levels = c(NMLL + 0.5, NMLL + 1.0, NMLL + 2.0, NMLL + 3.0),
   labels = c("1", "2", "4", "6"), labcex = 1.5,
    main = "Difference -2ll, alpha & beta, and =/- 2SE lines",
    xlab = "alpha",
	ylab = "beta")
	
# one can add lines at +/- 2 SEs to see what thsoe mean w.r.t the likelihood
lines(c(result$estimate[1]-2*SE[1], result$estimate[1]-2*SE[1]),
		c(y.beta[1], y.beta[100]), lty="dotted", col = "blue")
lines(c(result$estimate[1]+2*SE[1], result$estimate[1]+2*SE[1]),
		c(y.beta[1], y.beta[100]), lty="dotted", col = "blue")
lines(c(x.alpha[1], x.alpha[100]),
		c(result$estimate[2]-2*SE[2], result$estimate[2]-2*SE[2]), 
				lty="dotted", col = "blue")
lines(c(x.alpha[1], x.alpha[100]),
		c(result$estimate[2]+2*SE[2], result$estimate[2]+2*SE[2]), 
				lty="dotted", col = "blue")


# but back to fitting: here it is with analytical derivatives
# (in programming these, a good trick is to use the commented-out
# print statements in there (below) and be sure the gradient
# and hessian are the sensible, as well as the same as those
# from the derivative free routine (above) at the end
# then, these values will be used to be sure the C++ is correct 
ll.probitD <- function(params, x, N, r) {
	alpha <- params[1]
	beta <- params[2]
	Y <- alpha + (beta*x)
	P <- pnorm(Y)
	Q <- 1-P
	l <- (-1)*sum(r*log(P) + ((N-r)*log(Q)))
# gradient and hessian from B&J pp. 53-56
	h <- (1/sqrt(2*pi))*exp((-0.5)*Y*Y)
	p <- r/N
	dldalpha <- (-1)*sum((N*(p - P)*h)/(P*Q))
	dldbeta <- (-1)*sum((N*(p - P)*h*x)/(P*Q))
	d2ldalpha2 <- (-1)*(sum(((N*h*(p-P))/(P*Q))*
							((h/Q)-(h/P)-Y)) - sum((N*h*h)/(P*Q)))
	d2ldalphabeta <- (-1)*(sum(((N*h*(p-P))/(P*Q))*
							((h/Q)-(h/P)-Y)*x) - sum((N*h*h*x)/(P*Q)))
	d2ldbeta2 <- (-1)*(sum(((N*h*(p-P))/(P*Q))*
							((h/Q)-(h/P)-Y)*x*x) - sum((N*h*h*x*x)/(P*Q)))
	
	attr(l, "gradient") <- c(dldalpha, dldbeta)
	attr(l, "hessian") <- matrix(c(d2ldalpha2, d2ldalphabeta,
							d2ldalphabeta, d2ldbeta2), 2, 2)
#	print(params)
#	print(l)
	return(l)

 }

params <- c(0,5)

system.time( result <- nlm(f=ll.probitD,p=params,hessian=TRUE,x=x,N=N,r=r) )
print(result)
SE <- sqrt(diag(solve(result$hessian)))
print(SE)