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