### R code from vignette source '../confidence_intervals_and_hypothesis_testing_chapter/confidence_intervals_and_hypothesis_testing.Rnw'

###################################################
### code chunk number 1: preliminaries
###################################################
library(languageR)
library(bayesm)
my.cols <- c("magenta","green")
# first, function to fill area under curve
fill <- function(x1, x2, f,n=50, color="gray") {
  xvals <- seq(x1, x2, by=(x2-x1)/n)
  yvals <- f(xvals)
  #print(xvals)
  #print(yvals)
  #print(length(xvals))
  polygon(c(xvals, rev(xvals)), c(rep(0, n+1), rev(yvals)), col=color,border=NA)
  lines(c(xvals[1], xvals[length(xvals)]), c(0,0), col=color)
}
# second, function to prettify S scalar expression output
latex.scalar <- function(x) {
  # TODO
  x
}
options(lineWidth=55)


###################################################
### code chunk number 2: HPDinterval
###################################################
b1 <- 5
b2 <- 29
alpha <- 0.05
p <- seq(0,1,by=0.01)
old.par <- par(mar=c(5,4,4,5)+0.1)
plot(p, dbeta(p, b1, b2),type="l",ylab=expression(p(paste(pi, "|", bold(y)))),xlab=expression(pi),lwd=1,ylim=c(0,8))
minlength <- 1
p1 <- 0
p2 <- 1
for(p1.test in seq(0,1,by=0.01)) {
  mass.left <- pbeta(p1.test, b1, b2)
  if(mass.left > alpha)
    break
  p2.test <- qbeta(mass.left + 1 - alpha, b1, b2)
  if(p2.test - p1.test < minlength) {
    minlength <- p2.test - p1.test
    p1 <- p1.test
    p2 <- p2.test
  }
}
lines(c(p1,p1), c(0,dbeta(p1, b1, b2)))
lines(c(p2,p2), c(0,dbeta(p2, b1, b2)))
#fill(0,p1,function(x) dbeta(x, b1, b2))
fill(p1,p2,function(x) dbeta(x, b1, b2))


###################################################
### code chunk number 3: symmetricBayesianInterval
###################################################
b1 <- 5
b2 <- 29
alpha <- 0.05
p <- seq(0,1,by=0.01)
old.par <- par(mar=c(5,4,4,5)+0.1)
plot(p, dbeta(p, b1, b2),type="l",ylab=expression(p(paste(pi, "|", bold(y)))),xlab=expression(pi),lwd=1,ylim=c(0,8))
p1 <- qbeta(alpha/2, b1,b2)
p2 <- qbeta(1 - alpha/2, b1,b2)
lines(c(p1,p1), c(0,dbeta(p1, b1, b2)))
lines(c(p2,p2), c(0,dbeta(p2, b1, b2)))
#fill(0,p1,function(x) dbeta(x, b1, b2))
fill(p1,p2,function(x) dbeta(x, b1, b2))


###################################################
### code chunk number 4: confidence_intervals_and_hypothesis_testing.Rnw:375-378
###################################################
pyH1 <- choose(6,4) * 0.5^6
pyH2 <- 0.5 * choose(6,4) * (1/3)^2 * (2/3)^2 * ( (1/3)^2 + (2/3)^2)
py <- (pyH1 + pyH2) / 2


###################################################
### code chunk number 5: confidence_intervals_and_hypothesis_testing.Rnw:469-471
###################################################
pyH3 <- choose(6,4) * beta(5,3)
pH1y <- pyH1 / (pyH1 + pyH3)


###################################################
### code chunk number 6: probabilityDensitiesBandPVOT
###################################################
mu1 <- 0
mu2 <- 50
sigma <- 12
x <- seq(-20,80,by=0.1)
y1 <- dnorm(x,mu1,sigma)
y2 <- dnorm(x,mu2,sigma)
plot(y1 ~ x,type="l",col=my.cols[1],xlab="VOT",ylab="Probability density",lwd=2)
lines(y2 ~ x, col=my.cols[2],lwd=2,lty=2)
text(c(-10,60),c(0.031,0.031),c("/b/","/p/"),col=my.cols[1:2])
x1 <- 27.25
x2 <- 26.75
mean.x <- mean(c(x1,x2))
arrows(x1,0,x1,dnorm(x1,mu1,sigma)-0.00025,length=0.1,col=my.cols[1])
arrows(x2,0,x2,dnorm(x2,mu2,sigma)-0.00025,length=0.1,col=my.cols[2],lty=2)
points(mean.x,-0.00025,pch=19,cex=0.65)
axis(1,at=27,lwd=2,cex=0.5,tck=-0.01,font=6,cex.axis=0.75,padj=-1)


###################################################
### code chunk number 7: bayesPhonemeDiscriminationPosteriorProbs
###################################################
f <- function(x) 1/(1+exp( ((x-mu1)^2 - (x-mu2)^2 ) / (2*sigma^2)))
plot(f(x) ~ x, type="l",ylab="Posterior probability of /b/",xlab="VOT",lwd=2)
points(mean.x,f(mean.x),pch=19)


###################################################
### code chunk number 8: VOTVarianceManipulation
###################################################
sigma <- 8
y1 <- dnorm(x,mu1,sigma)
y2 <- dnorm(x,mu2,sigma)
plot(y1 ~ x,type="l",col="magenta",xlab="VOT",ylab="Probability density",lwd=1)
lines(y2 ~ x, col="green",lwd=1,lty=2)
sigma <- 14
y1 <- dnorm(x,mu1,sigma)
y2 <- dnorm(x,mu2,sigma)
lines(y1 ~ x,col="magenta",lwd=3)
lines(y2 ~ x, col="green",lwd=3,lty=2)
text(c(-10,60),c(0.035,0.035),c("[b]","[p]"),col=c("magenta","green"))


###################################################
### code chunk number 9: VOTVarianceIdealResponses
###################################################
sigma <- 14
y1 <- dnorm(x,mu1,sigma)
y2 <- dnorm(x,mu2,sigma)
discrim <- y1/(y1+y2)
plot(discrim ~ x, type="l",lwd=3,xlab="VOT",ylab="Posterior probability of /b/")
sigma <- 8
y1 <- dnorm(x,mu1,sigma)
y2 <- dnorm(x,mu2,sigma)
discrim <- y1/(y1+y2)
lines(discrim ~ x, lwd=1)


###################################################
### code chunk number 10: VOTVarianceEmpiricalResponses
###################################################
dat <- read.table("../data/clayards-etal-2008/cognition_trialbytrial_resp.txt",header=T)
rates <- with(subset(dat,VOT >= -20),aggregate(list(response=response),list(condition=condition,VOT=VOT),mean,na.rm=T))
plot(1-response ~ VOT, data=subset(rates,condition=="N"),type="b",col="black",lwd=1,pch='.',ylab="Proportion response /b/",xlim=c(-20,80))
lines(1-response ~ VOT, data=subset(rates,condition=="W"),type="b",col="black",lwd=3,pch=19,cex=0.5)


###################################################
### code chunk number 11: confidence_intervals_and_hypothesis_testing.Rnw:918-922
###################################################
pb <- read.table("../data/peterson_barney_data/pb.txt",header=T)
eh <- subset(pb,Vowel=="eh")
half.length <- round(sd(eh[["F1"]]) / sqrt(length(eh[["F1"]])) * (-1 * qt(0.025,length(eh[["F1"]]))),2)
eh.mean <- round(mean(eh[["F1"]]),2)


###################################################
### code chunk number 12: tBasedConfidenceInterval
###################################################
n <- 5
x <- seq(-4,4,by=0.01)
old.par <- par(mar=c(5,6,4,2)+0.1)
plot(x,dt(x,n),type="l",xlab=expression(frac(hat(mu) - mu, sqrt (S^2 / n))), ylab=expression(p(frac(hat(mu) - mu, sqrt (S^2 / n)))))
low.cutoff <- qt(0.025,n)
x.low <- subset(x,x < low.cutoff)
polygon(c(-4,x.low,low.cutoff),c(0,dt(x.low,n),0),col="lightgrey")
high.cutoff <- qt(0.975,n)
x.high <- subset(x,x >  high.cutoff)
polygon(c(high.cutoff,x.high,4),c(0,dt(x.high,n),0),col="lightgrey")
par(old.par)


###################################################
### code chunk number 13: epsilonF1
###################################################
#pb <- read.table("../data/peterson_barney_data/pb.txt",header=T)
#eh <- subset(pb,Vowel=="eh")
#length(eh[,1]) # 152 data points
#mean(eh$F1)
#sd(eh$F1)
hist(eh$F1,breaks=20,prob=T)
x <- seq(350,900,by=1)
lines(x,dnorm(x,590.7,97.1))


###################################################
### code chunk number 14: sigPowerTradeoff
###################################################
x <- seq(0,45,by=1)
plot(x,dbinom(x, 45, 0.5),type="l",ylim=c(0,0.2),xlab=expression(r),ylab=expression(P(r)),cex.lab=1.2,cex.axis=1.2)
legend(15, 0.2, c(expression(paste(pi[0],"=0.5")),expression(paste(pi[A],"=0.75")),expression(paste(pi[A],"=0.25"))), col=c("black","green","magenta"),lwd=3,cex=1.2)
fill(18-0.2,18+0.2,function(x) rep(dbinom(18,45,0.5),length(x)),color="darkgrey")
fill(27-0.2,29+0.2,function(x) dbinom(floor(x), 45, 0.5) + (dbinom(ceiling(x), 45, 0.5) - dbinom(floor(x), 45, 0.5)) * (x - floor(x)))
lines(x,dbinom(x, 45, 0.25),type="l",lwd=3,col="magenta")
lines(x,dbinom(x, 45, 0.75),type="l",lwd=3,col="green")
lines(x,dbinom(x, 45, 0.5),type="l",lwd=3,ylim=c(0,0.25),xlab=expression(r),ylab=expression(P(r)))


###################################################
### code chunk number 15: sigPowerTradeoffRatio
###################################################
old.par <- par(mar=c(5,5.5,4,2)+0.1)
plot(x,log(dbinom(x, 45, 0.5)/dbinom(x,45,0.75)),type="l",lwd=3,col="green",xlab=expression(r),
  ylab=expression(log(frac(paste("P(r | ", H[A], ")"),paste("P(r | ", H[0],"))")))),cex.axis=1.2,cex.lab=1.2) # FIX THE PROB SYMBOLS
lines(x,log(dbinom(x, 45, 0.5)/dbinom(x,45,0.25)),type="l",lwd=3,col="magenta")
legend(15, 30, c(expression(paste(pi[A],"=0.75")),expression(paste(pi[A],"=0.25"))), col=c("green","magenta"),lwd=3,cex=1.2)
par(old.par)


###################################################
### code chunk number 16: binomTestCorrectRejectionRegion
###################################################
x <- seq(0,45,by=1)
plot(x,dbinom(x, 45, 0.5),type="l",ylim=c(0,0.15),xlab=expression(r),ylab=expression(P(r)),lwd=3,lab=c(10,5,7),cex.axis=1.2,cex.lab=1.2)
fill(0,15+0.2,function(x) dbinom(floor(x), 45, 0.5) + (dbinom(ceiling(x), 45, 0.5) - dbinom(floor(x), 45, 0.5)) * (x - floor(x)))
fill(30-0.2,45,function(x) dbinom(floor(x), 45, 0.5) + (dbinom(ceiling(x), 45, 0.5) - dbinom(floor(x), 45, 0.5)) * (x - floor(x)))
fill(15,16,function(x) dbinom(floor(x), 45, 0.5) + (dbinom(ceiling(x), 45, 0.5) - dbinom(floor(x), 45, 0.5)) * (x - floor(x)),col="black")