## A psychophysical experiment - 
## This runs 330 trials, i.e., requires 330 responses from observer
runSampleExperiment(DisplayTrial = "DisplayOneTrial", 
  DefineStimuli = "DefineMyScale")

## KK's results - from running the experiment above 3 times on one observer
data("kk1")
data("kk2")
data("kk3")
kk <- SwapOrder(rbind(kk1, kk2, kk3))

## MLDS - using glm (default)
kk.mlds <- mlds(kk)
summary(kk.mlds)

##  MLDS - using optim
kkopt.mlds <- mlds(kk, method = "optim",
  opt.init = c(seq(0, 1, length.out = 11), 0.2))
summary(kkopt.mlds)

##  MLDS - using alternative links
kk.mlds.logit <- mlds(kk, lnk = "logit")
kk.mlds.cauchit <- mlds(kk, lnk = "cauchit")

##  estimating deviance accounted for
(kk.mlds$obj$null.deviance - deviance(kk.mlds$obj))/kk.mlds$obj$null.deviance

## bootstrap estimate of empirical cdf of residuals
kk.diag.prob <- binom.diagnostics(kk.mlds, 10000)
str(kk.diag.prob)
kk[residuals(kk.mlds$obj) < -2.5, ]  ##which trials correspond to outlier residuals?

## Figure 1
## The bivariate point distributions were generated with the
## mvrnorm() function from the package \pkg{MASS}.
library("MLDS")
rc <- DefineMyScale()
N <- 100
ptSize <- 1
opar <- par(mfrow = c(3, 4))
for (ix in 1:11) {
	covm <- matrix(c(1, rep(rc[ix], 2), 1), 2, 2)
	xy <- MASS:::mvrnorm(n = N, rep(0, 2), covm, 
		empirical = TRUE)
	plot(xy, axes = FALSE, xlab = "", ylab = "", pty = "s",
		cex = ptSize, pch = 16, col = "black",
		xlim = c(-4, 4), ylim = c(-4, 4),
		main = substitute(paste(r, " = ", rc), 
	list(rc = round(rc[ix], 2))), cex.main = 1.75  ) 
	}
data("kk3")
par(mar = c(7, 5, 4, 2) + 0.1)
plot(mlds(SwapOrder(kk3), method = "optim", 
		opt.init = c(seq(0, 1, len = 11), 0.2)), type = "b",
		xlab = "r", ylab = "Difference Scale Value",
		cex.axis = 1.1, cex = 1.25, cex.lab = 1.4, 
		las = 1)	
par(opar)

## Figure 2 - a sample trial from the experiment
opar <- par(mfrow = c(2, 2))
DisplayOneTrial(c(0.6, 0.9, 0.1, 0.3))
par(opar)

## Figure 3 - the 6-point test, schematically
opar <- par(xaxt = "n", yaxt = "n", bty = "n")
plot(c(0, 14), c(0, 1), type = "n", xlab = "", ylab = "")
lines(c(1, 12), c(0.5, 0.5), lwd = 5)
points(c(3.25, 5.5, 7, 9.15, 10.25, 11.25), rep(0.5, 6),
	pch = 16, cex = 1.55)
text(3.25, 0.2, "a", cex = 3)
text(5.5, 0.2, "b", cex = 3)
text(7, 0.2, "c", cex = 3)
text(9.15, 0.2, "a'", cex = 3)
text(10.25, 0.2, "b'", cex = 3)
text(11.25, 0.2, "c'", cex = 3)
par(opar)

## Figure 4
### 6-point condition on circle
opar <- par(lwd = 3, pty = "s")
th <- seq(0, pi, len = 200)
plot(cos(th), sin(th), type = "l", axes = FALSE,
	xlab = "", ylab = "", xlim = c(-1.3, 1.3), ylim = c(-1.3, 1.3))
th <- seq(pi, 1.2*pi/2, len = 3)
lines(cos(th), sin(th))
rd <- c(pi,1.2*pi/2)
lines(cos(rd), sin(rd), col = "red", lty = 2)
points(cos(th), sin(th), pch = 16, col = "black")
text(1.2 * cos(th), 1.2 * sin(th), c("a", "b", "c"), cex = 3)
th <- seq(pi/2,0.2*pi/2, len = 3)
points(cos(th), sin(th), pch = 16, col = "black")
lines(cos(th), sin(th))
rd <- c(pi/2,0.2*pi/2)
lines(cos(rd), sin(rd), col = "red", lty = 2)
text(1.2 * cos(th), 1.2 * sin(th), c("a'", "b'", "c'"), cex = 3)
mtext("a", line = 1, at = -1, adj = 0, cex = 4)
par(opar)
### 6-point condition on ellipse
opar <- par(lwd = 3, pty = "s")
asr <- 0.475
dd <- function(th, x, y, d) {
	sqrt((cos(th) - x)^2 + (asr *  sin(th) - y)^2)  - d }

th <- seq(0, pi, len = 200)
plot(cos(th), asr * sin(th), type = "l", axes = FALSE,
	xlab = "", ylab = "", xlim = c(-1.3, 1.3), 
	ylim = c(-1.3, 1.3))
th <- seq(pi, 5 * pi/8, len = 3) #- 1  #points a, b, c
dst <- sqrt(rowSums((apply(cbind(cos(th),  asr * sin(th)), 2, diff))^2))
points(cos(th),  asr * sin(th), pch = 16, col = "black")
lines( cos(th), asr * sin(th))
lines(cos(th[-2]), asr *  sin(th[-2]), col = "red", lty = 2)
text(1.225 * cos(th), 1.225 *  asr * sin(th), c("a", "b", "c"), cex = 3)
th2 <- th[3] - pi/16  # the point of this is to equate the chords
th2 <- c(th2, 
	uniroot(dd, c(th2, th2 - pi/6), 
	x =  cos(th2), y = asr * sin(th2), d = dst[1])$root
	)
th2 <- c(th2, 
	uniroot(dd, c(th2[2], th2[2] - pi^2/10), 
		x = cos(th2[2]), y =  asr * sin(th2[2]), d = dst[2])$root
		)
points(cos(th2),  asr * sin(th2), pch = 16, col = "black")
lines(cos(th2),  asr * sin(th2))
lines(cos(th2[-2]),  asr * sin(th2[-2]), col = "red", lty = 2)
text(1.25 * cos(th2), 1.25 *  asr * sin(th2), c("a'", "b'", "c'"), cex = 3)
ac0 <- sum(rowSums(diff(cbind(cos(th[-2]), asr * 	sin(th[-2])))^2))
ac1 <- sum(rowSums(diff(cbind(cos(th2[-2]), asr * 
	sin(th2[-2])))^2))
txt <- paste("ac = ", round(ac0, 2), ", a'c' = ", 
	round(ac1, 2), sep = "")
text(0, -0.75, txt, cex = 3)
mtext("b", line = 1, at = -1, adj = 0, cex = 4)
par(opar)

## Figure 5
## The standard deviations of the bootstrap replications were calculated
## in advance for efficiency.
kk.sd <- c(0, 0.0245, 0.0252, 0.02783, 0.02461, 0.02039, 
	0.01913, 0.01859, 0.01785, 0.01862, 0)
opar <- par(mfrow = c(1, 2))
plot(kk.mlds, standard.scale = TRUE, xlab = expression(r), 
	ylab = "Difference Scale Value", 
	cex.lab = 1.5, cex.axis = 1.5, las = 1)
psc <- with(kk.mlds, pscale/pscale[length(pscale)])
segments(kk.mlds$stimulus, psc + 2 * kk.sd, 
	kk.mlds$stimulus, psc -2 * kk.sd)
lines(kkopt.mlds$stimulus, kkopt.mlds$pscale)
plot(kk.mlds$stimulus^2, kk.mlds$pscale/max(kk.mlds$pscale), 
	xlab = expression(r^2), ylab = "Difference Scale Value", 
	cex.lab = 1.5, cex.axis = 1.5, las = 1, xlim = c(0, 1))
abline(0, 1, lty = 2, lwd = 2)
lines(kkopt.mlds$stimulus^2, kkopt.mlds$pscale)
par(opar)

## Figure 6
dd <- data.frame(resp = as.factor(mlds(SwapOrder(kk2))$obj$data$resp), 
		fitted = fitted(mlds(SwapOrder(kk2))), dataset = "kk2")
ee <- data.frame(resp = as.factor(x.mlds$obj$data$resp), 
		fitted = fitted(x.mlds), dataset = "AutumnLab")
dd <- rbind(dd, ee)
library("ggplot")
p <- ggplot(dd, . ~ dataset, aes = list(x = resp, y = fitted))
p$xlabel <- "Response"
p$ylabel <- "Fitted Probability"
strnames <- c(kk2 = "kk2", Autumn = "AutumnLab")
ggopt(strip.text = function(variable, value) 
	strnames[value]     )
print(ggpoint(p, size = 3, colour = rgb(0, 0, 0, 0.25)))

## Figure 7 - compare the empirical cdf's of residuals 
##   and the number of runs
kk.mlds.logit <- mlds(kk, lnk = "logit")
kk.mlds.cauchit <- mlds(kk, lnk = "cauchit")
###each of the next 3 lines will take about 15 min.
kk.diag.prob <- binom.diagnostics(kk.mlds, 10000)
kk.diag.log <- binom.diagnostics(kk.mlds.logit, 10000)
kk.diag.cau <- binom.diagnostics(kk.mlds.cauchit, 10000)  
plot(kk.diag.prob, cex = 0.5,  main = "probit link", breaks = 20)
plot(kk.diag.log, cex = 0.5, main = "logit link", breaks = 20)
plot(kk.diag.cau, cex = 0.5, main = "cauchit link", breaks = 20)

## Figure 8 - 6-point test
### expect to wait about 15 minutes and some warnings
kk.6pt <- simu.6pt(kk.mlds, 10000)   
str(kk.6pt)
hist(kk.bt$boot.samp, xlab = "log(Likelihood)", cex.lab = 1.5)
abline(v = kk.bt$lik6pt, lwd = 3)