##### script for slidedeck cA
##### 
##### Slide deck for Approximating the Binomial
##### 



### Preamble

ito1 = "#FAFAFF"  ## background
ito2 = "#00AAE4"  ## selected






#####
###  Graphic: Binom(5,0.500) and N(2.5, 1.25)


n=5; p=0.50
xx = seq(-5,10, length=1e4)
yb = pbinom(xx, size=n, prob=p)
yn = pnorm(xx, m=n*p,s=sqrt(n*p*(1-p)))



# Binomial pmf
par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,1,1)+0.75)
par(font.lab=2,cex.lab=1.1)
par(xpd=NA)

plot.new()
plot.window( xlim=c(-2,7), ylim=c(0,1.01))

axis(1, at=c(-5,0:5,10), label=c(NA,0:5,NA))
axis(2)
title(xlab="Random Variable", line=2.5)
title(ylab="Probability", line=2.75)

lines(0:5, dbinom(0:5, size=n, prob=p), col="red", type="h", lwd=2)
points(0:5, dbinom(0:5, size=n, prob=p), col="red", pch=16)



# Binomial CDF
par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,1,1)+0.75)
par(font.lab=2,cex.lab=1.1)
par(xpd=NA)

plot.new()
plot.window( xlim=c(-2,7), ylim=c(0,1))

points(xx,yb, col="red", pch=".", cex=3)

x2 = seq(0,5)
y2 = pbinom(x2-0.01,n,p)
points(x2,y2, col="red", pch=21, bg="white")

x2 = seq(0,5)
y2 = pbinom(x2,n,p)
points(x2,y2, col="red", pch=16)

axis(1, at=-2:7)
axis(2)
title(xlab="Random Variable", line=2.5)
title(ylab="Cumulative Probability", line=2.75)





# Normal CDF

par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,1,1)+0.75)
par(font.lab=2,cex.lab=1.1)
par(xpd=NA)

plot.new()
plot.window( xlim=c(-2,7), ylim=c(0,1))

axis(1, at=-2:7)
axis(2)
title(xlab="Random Variable", line=2.5)
title(ylab="Cumulative Probability", line=2.75)

lines(xx,yn, col="blue", lwd=2)






# Both, together

par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,1,1)+0.75)
par(font.lab=2,cex.lab=1.1)
par(xpd=NA)

plot.new()
plot.window( xlim=c(-2,7), ylim=c(0,1))

axis(1, at=-2:7)
axis(2)
title(xlab="Random Variable", line=2.5)
title(ylab="Cumulative Probability", line=2.75)

points(xx,yb, col="red", pch=".", cex=3)
lines(xx,yn, col="blue", lwd=2)

x2 = seq(0,5)
y2 = pbinom(x2-0.01,n,p)
points(x2,y2, col="red", pch=21, bg="white")

x2 = seq(0,5)
y2 = pbinom(x2,n,p)
points(x2,y2, col="red", pch=16)










#####
###  The Examples


### Ex1: Cheating
n = 50
p = 1/50

pnorm(1.5, m=n*p, sqrt(n*p*(1-p)) )
pbinom(1, size=n, prob=p)




### Ex2: Births
n = 100
p = 0.500

1-pnorm(55.5, n*p, sqrt(n*p*(1-p)) )
1-pbinom(55,n,p)





### Ex3: Stats Examination

n=50
p=0.10

pn = pnorm(4.5, m=n*p, s=sqrt(n*p*(1-p)))
pb = pbinom(4, n,p)
pn;pb

pn-pb
(pn-pb)/pb * 100

pbinom(3, n,p)



### Ex4: Toothpaste

n=400
p=0.750

pn = pnorm(300.5, m=n*p, s=sqrt(n*p*(1-p))) - 
	 pnorm(299.5, m=n*p, s=sqrt(n*p*(1-p)))

pb = dbinom(300, n,p)

pn;pb

pn-pb
(pn-pb)/pb * 100

pbinom(3, n,p)

















#####
###  Continuity Correction Examples

### Ex1: Cheating Graphic

mm=1; ss=sqrt(0.98)
rptChg = 0.2

xx = seq(-10,10, length=1e4)
yy = dnorm(xx, mm,ss)




par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,1,1)+0.75)
par(font.lab=2,cex.lab=1.1)

plot.new()
plot.window( xlim=c(-0.9,4.9), ylim=c(0,0.42))

axis(1, at=c(-10,-1:4,100), label=c(NA,-1:4,NA))
axis(2)
title(xlab="Random Variable (Cheaters)", line=2.5)

polygon( x=c(xx,rev(xx)), y=c(yy, rep(0,length(yy))), col=ito1 )




# Leftmost
thisX = seq(-5,-1.5,length=1e4)
thisY = dnorm(thisX, mm,ss)
polygon( x=c(thisX,rev(thisX)), y=c(thisY, rep(0,length(thisY))), col=rgb(0,0,1,rptChg) )

# -0.5
thisX = seq(-5,-0.5,length=1e4)
thisY = dnorm(thisX, mm,ss)
polygon( x=c(thisX,rev(thisX)), y=c(thisY, rep(0,length(thisY))), col=rgb(0,0,1, rptChg) )


# 0.5
thisX = seq(-5,0.5,length=1e4)
thisY = dnorm(thisX, mm,ss)
polygon( x=c(thisX,rev(thisX)), y=c(thisY, rep(0,length(thisY))), col=rgb(0,0,1, rptChg) )


# 1.5
thisX = seq(-5,1.5,length=1e4)
thisY = dnorm(thisX, mm,ss)
polygon( x=c(thisX,rev(thisX)), y=c(thisY, rep(0,length(thisY))), col=rgb(0,0,1, rptChg) )


# 2.5
thisX = seq(-5,2.5,length=1e4)
thisY = dnorm(thisX, mm,ss)
polygon( x=c(thisX,rev(thisX)), y=c(thisY, rep(0,length(thisY))), col=rgb(0,0,1, rptChg) )


# 3.5
thisX = seq(-5,3.5,length=1e4)
thisY = dnorm(thisX, mm,ss)
polygon( x=c(thisX,rev(thisX)), y=c(thisY, rep(0,length(thisY))), col=rgb(0,0,1, rptChg) )


# 4.5
thisX = seq(-5,4.5,length=1e4)
thisY = dnorm(thisX, mm,ss)
polygon( x=c(thisX,rev(thisX)), y=c(thisY, rep(0,length(thisY))), col=rgb(0,0,1, rptChg) )



axis(1, at=seq(-10.5,10.5,1), label=NA)
lines(xx,yy, lwd=2, col="black")











### Ex2: Births Graphic


rptChg = 0.2

n=100; p=0.500
mm=n*p; ss=n*p*(1-p)

xx = seq(0,100, length=1e4)
yy = dnorm(xx, mm, sqrt(ss) )

xxxx = seq(55,100, length=1e4)
yyyy = dnorm(xxxx, mm, sqrt(ss) )




par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,2,1,1)+0.75)
par(font.lab=2,cex.lab=1.1)
par(xpd=NA)

plot.new()
plot.window( xlim=c(25,75), ylim=c(0,0.08))

axis(1, at=c(-10,seq(0,200,10),300), label=c(NA,seq(0,200,10),NA))
axis(1, at=seq(25,75,5), label=NA)
axis(2, at=seq(0,2,0.02))
title(xlab="Random Variable (Births of Girls)", line=2.5)

polygon( x=c(xx,rev(xx)), y=c(yy, rep(0,length(yy))), col=ito1 )
polygon( x=c(xxxx,rev(xxxx)), y=c(yyyy, rep(0,length(yyyy))), col=ito2 )

lines(xx,yy, lwd=2, col="black")

x3 = seq(0,100,5)
y3 = dbinom(x3, n, p)
points(x3,y3, pch=21, bg="white")





















#####
###  The Animation


# Animation of the previous
library(animation)

ani.options(interval=0.25)

# Begin animation loop
saveGIF({
  # Adjust things a little
  
  par(bg="transparent")
  par(xaxs="i", yaxs="i", xpd=TRUE)
  par(family="serif",las=1)
  par(mar=c(3,3,0,0)+0.75)
  par(font.lab=2,cex.lab=1.1)

for (n in 10:200) {
    

p=0.50
xx = seq(-100,200, length=1e4)
yb = pbinom(xx, size=n, prob=p)
yn = pnorm(xx, m=n*p,s=sqrt(n*p*(1-p)))



plot.new()
plot.window( xlim=c(-2,30), ylim=c(-0.1,0.1))

#axis(1, at=0:20)
abline(h=0)
axis(2, at=c(0.1,0.05,0.03,0.02)*c(-1,1))
title(xlab="Random Variable", line=2.5)
title(ylab="Difference in Cumulative Probabilities", line=2.75)

lines(xx-n*p+15,yn-yb, col="blue", lwd=2) ## Difference
text(x=0,y=0.1, label=paste("n =",n))
  }
},  movie.name = "cA-animation.gif", height=ht, width=wt)










# ===== ----- ===== -----
#
# This loop creates the many, many, many images used for 
# the animation in LaTeX.
#

saveLatex({
  # Adjust the margins a little

  par(xaxs="i", yaxs="i")
  par(family="serif", las=1)
  par(mar=c(3,3,0,3)+0.5)
  par(font.lab=2, cex.lab=1.1, cex.axis=0.8)
  par(cex=1.3)


for (n in 1:500) {

p=0.50
xx = seq(-100,1100, length=1e4)
yb = pbinom(xx, size=n, prob=p)
yn = pnorm(xx, m=n*p,s=sqrt(n*p*(1-p)))



plot.new()
plot.window( xlim=c(-2,30), ylim=c(-0.1,0.1))

theGrid = c(0.10,0.07,0.05,0.02,-0.02,-0.05,-0.07,-0.1)
axis(2, at=theGrid, label=paste(theGrid*100,"pp",sep="") )
axis(4, at=theGrid, label=paste(theGrid*100,"pp",sep="") )

abline(h=theGrid, lty=3)
abline(h=0, lwd=2, col="black")

title(xlab="The Random Variable", line=0.5)
mtext(side=2, text="Difference in Probabilities", line=2.5, cex=1.2, font=2, las=0)
mtext(side=4, text="Difference in Probabilities", line=2.5, cex=1.2, font=2, las=0)

lines(xx-n*p+15,yn-yb, col="blue", lwd=2) ## Difference
text(x=0,y=0.09, label=paste("n =",n))
  }
}, 	documentclass=NULL, latex.filename="cA-binAnim8", 
	ani.height=ht*1.1, ani.width=wt, img.name="cA-binAnim8_"

) ## End of LaTeX images










##### End of File