### R code from vignette source '~/WindowsC/HOME/rmh/hh.e2/hh2/iinf.tex'

###################################################
### code chunk number 1: iinf.tex:118-143
###################################################
hhpdf("Hyz6.pdf", height=8.5, width=7)

NTfun <- function(...)
   NTplot(..., main=NULL, xlab=NULL, ylab=NULL, prob.labels=FALSE,
          xhalf.multiplier=1.6, yhalf.multiplier=2, cex.top.axis=.7,
          scales=list(cex=.5, y=list(at=1)),
          xlim=c(-3, 4))

Nhr  <- NTfun(mean0=0, xbar= 1.8, alpha.right=.05,  alpha.left=0)
Nhl  <- NTfun(mean0=0, xbar=-1.8, alpha.right=0,    alpha.left=.05)
Nhlr <- NTfun(mean0=0, xbar= 1.8, alpha.right=.025, alpha.left=.025)

Ncl  <- NTfun(mean0=0, xbar= 1.8, alpha.right=0,    alpha.left=.05,  type="confidence")
Ncr  <- NTfun(mean0=0, xbar=-1.8, alpha.right=.05,  alpha.left=0,    type="confidence")
Nclr <- NTfun(mean0=0, xbar= 1.8, alpha.right=.025, alpha.left=.025, type="confidence")

print(Nhr,  tablesOnPlot=FALSE, position=c(.00, .630, .55, 1.000), more=TRUE)
print(Nhl,  tablesOnPlot=FALSE, position=c(.00, .315, .55,  .685), more=TRUE)
print(Nhlr, tablesOnPlot=FALSE, position=c(.00, .000, .55,  .370), more=TRUE)

print(Ncl,  tablesOnPlot=FALSE, position=c(.45, .630, 1.0, 1.000), more=TRUE)
print(Ncr,  tablesOnPlot=FALSE, position=c(.45, .315, 1.0,  .685), more=TRUE)
print(Nclr, tablesOnPlot=FALSE, position=c(.45, .000, 1.0,  .370), more=FALSE)

hhdev.off()


###################################################
### code chunk number 2: iinf.tex:369-394
###################################################
hhpdf("Hyp6.pdf", height=8.5, width=7)

NTfun <- function(...)
   NTplot(..., main=NULL, xlab=NULL, ylab=NULL, prob.labels=FALSE,
          xhalf.multiplier=1.6, yhalf.multiplier=2, cex.top.axis=.7,
          scales=list(cex=.5),
          xlim=c(0, 1), ylim=c(0, 4), distribution.name="binomial")

Nphr  <- NTfun(p0=.5, p.hat= .75, n=20, alpha.right=.05,  alpha.left=0)
Nphl  <- NTfun(p0=.5, p.hat= .25, n=20, alpha.right=0,    alpha.left=.05)
Nphlr <- NTfun(p0=.5, p.hat= .75, n=20, alpha.right=.025, alpha.left=.025)

Npcl  <- NTfun(p0=.5, p.hat= .75, n=20, alpha.right=0,    alpha.left=.05,  type="confidence")
Npcr  <- NTfun(p0=.5, p.hat= .25, n=20, alpha.right=.05,  alpha.left=0,    type="confidence")
Npclr <- NTfun(p0=.5, p.hat= .75, n=20, alpha.right=.025, alpha.left=.025, type="confidence")

print(Nphr,  tablesOnPlot=FALSE, position=c(.00, .630, .55, 1.000), more=TRUE)
print(Nphl,  tablesOnPlot=FALSE, position=c(.00, .315, .55,  .685), more=TRUE)
print(Nphlr, tablesOnPlot=FALSE, position=c(.00, .000, .55,  .370), more=TRUE)

print(Npcl,  tablesOnPlot=FALSE, position=c(.45, .630, 1.0, 1.000), more=TRUE)
print(Npcr,  tablesOnPlot=FALSE, position=c(.45, .315, 1.0,  .685), more=TRUE)
print(Npclr, tablesOnPlot=FALSE, position=c(.45, .000, 1.0,  .370), more=FALSE)

hhdev.off()


###################################################
### code chunk number 3: iinf.tex:497-501
###################################################
data(vocab)
hhcapture("vocab-stem.Rout", '
stem(vocab$score, scale=2)
')


###################################################
### code chunk number 4: iinf.tex:515-522
###################################################
hhcapture("vocab-t.Rout", '
vocab.t <- t.test(vocab$score, mu=10)
vocab.t
')
hhpdf("vocab.pdf", height=6, width=7)
NTplot(vocab.t, type="confidence", cex.prob=1)
hhdev.off()


###################################################
### code chunk number 5: iinf.tex:924-996
###################################################
data(teachers)
teachers$"English-Greek" <- teachers$English - teachers$Greek
teachers$sentence <- ordered(1:32, order(teachers$"English-Greek"))
teachers$smaller <- factor(sign(teachers$"English-Greek"),
                           labels=c("English\nsmaller error count",
                           "Greek\nsmaller error count"))
teachers[c(1,2,32),]

teachers.melt <- melt(cbind(teachers, sentence=1:32),
                      id.vars=c("sentence","English-Greek", "smaller"),
                      variable.name="language",
                      value.name="errors")
teachers.melt[c(1,2,32,33,63,64),]

hhpdf("teachers-dot.pdf", height=6.9, width=9)
tpg <- trellis.par.get("superpose.symbol")
tpg$pch[1:3] <- c(19,15,17)

AAAA <-
resizePanels(h=table(teachers$smaller),
dotplot(rep(teachers$sentence, 2) ~ errors | smaller,
        groups=language, data=teachers.melt,
        scales=list(y=list(relation="free")), pch=tpg$pch,
        par.strip.text=list(cex=1.1, lines=2),
        ylab="sentence",
        layout=c(1,2), between=list(y=1),
        key=list(text=
                   list(c("English for\nGreek Speakers",
                          "Greek for\nEnglish Speakers")),
          points=Rows(tpg, 1:2), padding.text=3,
          space="right",  border=TRUE)
        ))
## AAAA

BBBB <-
resizePanels(h=table(teachers$smaller),
             dotplot(sentence ~ `English-Greek` | smaller,
                     data=teachers, ref=0, pch=tpg$pch[3], col="gray30",
                     scales=list(y=list(relation="free")),
                     par.strip.text=list(cex=1.1, lines=2),
                     main="English - Greek",
                     layout=c(1,2), between=list(y=1)
                     )
             )
## BBBB + layer(panel.abline(v=0, col="gray60"))

CCCC <-
resizePanels(h=table(teachers$smaller),
             dotplot(sentence ~ sqrt(`English-Greek` + 17) | smaller,
                     data=teachers, ref=0, pch=tpg$pch[3], col="gray30",
                     par.strip.text=list(cex=1.1, lines=2),
                     scales=list(y=list(relation="free")),
                     main=expression(sqrt("English - Greek" + 17)),
                     layout=c(1,2), between=list(y=1))
                     )
## CCCC + layer(panel.abline(v=sqrt(17), col="gray60"))

AAAABBBBCCCC <-
resizePanels(h=table(teachers$smaller), w=c(30, 25, 25),
update(cbind(errors=AAAA, "English-Greek"=BBBB, "sqrt(17 +\nEnglish-Greek)"=CCCC),
       xlab=NULL, ylab=list("sentence", cex=1.2),
       par.strip.text=list(cex=1.1, lines=1.2),
       scales=list(
         y=list(relation="free"),
         x=list(limits=list(c(9, 42), c(-20, 20), c(0,6))))))
## AAAABBBBCCCC + layer(panel.abline(v=c(0,0,sqrt(17))[current.column()], col="gray60"))

xlab.top <- c(AAAABBBBCCCC$condlevels[[1]][1:2], expression(sqrt("English-Greek + 17")))
AAAABBBBCCCC$condlevels[[1]] <- c(" "," "," ")
AAAABBBBCCCC + layer(panel.abline(v=c(0,0,sqrt(17))[current.column()], col="gray60"))
grid.text(xlab.top, x=c(.25, .475, .695), y=.92)
hhdev.off()


###################################################
### code chunk number 6: iinf.tex:1043-1051
###################################################
hhcapture("teachersLeft.Rout", '
stem(teachers$"English-Greek")
t.test(teachers$"English-Greek")
')
hhcapture("teachersRight.Rout", '
stem(sqrt(teachers$"English-Greek" + 17), scale=.5)
t.test(sqrt(teachers$"English-Greek" + 17), mu=sqrt(17))
')


###################################################
### code chunk number 7: iinf.tex:1159-1184
###################################################
hhpdf("normalOneTailFig.pdf", height=7.2, width=9)
mean0 <- 1
delta <- 2
mean1 <- mean0 + delta
sd <- 3
alpha <- .05
beta <- .20
power <- 1-beta
n <- sd^2 * (qnorm(1-alpha) + qnorm(1-beta))^2 / delta^2
n
tmp <-
NTplot(mean0=mean0, mean1=mean1, sd=sd, n=n,
        key.axis.padding=6, cex.main=1.5,
        cex.top.axis=1.8, cex.prob=1.5,
        prob.labels=FALSE,
        digits.axis=5, digits.float=3, digits.left=3,
        xhalf.multiplier=.8, yhalf.multiplier=.6,
        zaxis=TRUE, z1axis=TRUE, cex.z=1)
print(tmp +
        layer({panel.segments(1, -.080, 3.000, -.080, col="gray60",    lwd=6) ## xbar scale
               panel.segments(1, -.098, 2.323, -.098, col="lightblue", lwd=6) ## z scale
               panel.segments(3, -.122, 2.323, -.122, col="#FF808B",   lwd=6) ## z1 scale
              }, under=TRUE),
      cex.table=1.2, digits=5)
hhdev.off()


###################################################
### code chunk number 8: iinf.tex:1256-1269
###################################################
hhcapture("sampleSizeA.Rout", '
## one sided
alpha <- .05
power <- .80
beta <- 1-power
delta <- 1
sd <- 2

## Approximation using formula assuming normal is appropriate
sd^2*(qnorm(1-alpha) + qnorm(1-beta))^2 / delta^2
## [1] 24.73
## n is slightly smaller with the normal assumption.
')


###################################################
### code chunk number 9: iinf.tex:1285-1294
###################################################
hhcapture("sampleSizeB.Rout", '
## solve using power.t.test
PTT <-
power.t.test(delta=delta, sd=sd, sig.level=alpha, power=power,
             type="one.sample", alternative="one.sided")
PTT
NTplot(PTT, zaxis=TRUE)  ## static plot
## NTplot(PTT, zaxis=TRUE), shiny=TRUE)  ## dynamic plot
')


###################################################
### code chunk number 10: iinf.tex:1310-1336
###################################################
hhcapture("sampleSizeC.Rout", '
## solve manually with t distribution.  Use ncp for alternative.
n0 <- 30 ## pick an n0 for starting value
t.critical <- qt(1-alpha, df=n0-1)
t.critical
## [1] 1.699

## a series of n values
nn <- 23:30
names(nn) <- nn
nn

## find the power for a series of n values for the specified critical value
pt(t.critical, df=nn-1, ncp=delta/(sd/sqrt(nn)), lower=FALSE)
##     23     24     25     26     27     28     29     30
## 0.7568 0.7722 0.7868 0.8006 0.8136 0.8258 0.8374 0.8483

## recalculate critical value with new n=26
t.critical <- qt(1-alpha, df=26-1)
t.critical
## find the power for a series of n values for the new critical value
pt(t.critical, df=nn-1, ncp=delta/(sd/sqrt(nn)), lower=FALSE)
##     23     24     25     26     27     28     29     30
## 0.7540 0.7695 0.7842 0.7981 0.8112 0.8235 0.8352 0.8461
## conclude n between 26 and 27
')


###################################################
### code chunk number 11: iinf.tex:1394-1420
###################################################
tmp26 <- NTplot(NTplot(PTT), main=NA,
                key.axis.padding=6, cex.main=1.5,
                cex.top.axis=1.8, cex.prob=1.5,
                prob.labels=FALSE,
                digits.axis=5, digits.float=3, digits.left=3,
                xhalf.multiplier=.8, yhalf.multiplier=.6,
                zaxis=TRUE, cex.z=1)
## tmp26

tmp30 <- NTplot(NTplot(PTT), n=30, df=29, main=NA,
                key.axis.padding=6, cex.main=1.5,
                cex.top.axis=1.8, cex.prob=1.5,
                prob.labels=FALSE,
                digits.axis=5, digits.float=3, digits.left=3,
                xhalf.multiplier=.8, yhalf.multiplier=.6,
                zaxis=TRUE, cex.z=1)
## tmp30

hhpdf("sampleSizeFig26.pdf", height=8, width=6) ## density, power, beta by sample size
print(update(tmp26, ylab=NULL), cex.table=1.2)
hhdev.off()

hhpdf("sampleSizeFig30.pdf", height=8, width=6) ## density, power, beta by sample size
print(update(tmp30, ylab=NULL), cex.table=1.2)
hhdev.off()



###################################################
### code chunk number 12: iinf.tex:1586-1592
###################################################
hhcapture("fairDice.Rout", '
dice <- sample(rep(1:6, c(3,7,5,8,1,6)))
dice
table(dice)
chisq.test(table(dice))
')


###################################################
### code chunk number 13: iinf.tex:1645-1665
###################################################
hhcapture("ex.chisq.Rout", '
Observed <- c(13, 18, 20, 18, 6, 5)
names(Observed) <- 0:5
## binomial proportion p=.4 is specified
Expected <- dbinom(0:5, size=5, p=.4)*80
names(Expected) <- 0:5
chisq.test(Observed, p=Expected, rescale.p=TRUE)


## binomial proportion p is calculated from the observations
p <- sum(Observed * (0:5)/5)/sum(Observed)
p
Expected <- dbinom(0:5, size=5, p=p)*80
names(Expected) <- 0:5
WrongDF <- chisq.test(Observed, p=Expected, rescale.p=TRUE)
WrongDF
c(WrongDF$statistic, WrongDF$parameter)
## correct df and p-value
pchisq(WrongDF$statistic, df=WrongDF$parameter - 1, lower=FALSE)
')


###################################################
### code chunk number 14: iinf.tex:1816-1836
###################################################
hhpdf("iinf-f-qqnorm.pdf", height=7.2, width=9)

qqnorm.dat <- data.frame(
   normal=rnorm(100,0,1),
   uniform=runif(100,-3,3))
qqnorm.dat <- within(qqnorm.dat, {
  `heavy-tailed` <- normal/uniform
  `thin-tailed` <- normal*exp(-abs(normal))
  `positively skewed` <- normal^2
  `negatively skewed` <- -`positively skewed`
  quantiles <- qnorm(ppoints(length(normal)))
})
qqnorm.dat <- data.frame(lapply(qqnorm.dat, sort), check.names=FALSE) ## sort on quantile is redundant
head(qqnorm.dat)

xyplot(normal + uniform + `negatively skewed` + `positively skewed` +
      `thin-tailed` + `heavy-tailed` ~ quantiles, data=qqnorm.dat, outer=TRUE,
      ylab="Observed values", xlab="Quantiles of Standard Normal",
      scales=list(y=list(relation="free")), between=list(x=1, y=1), layout=c(2,3))
hhdev.off()


###################################################
### code chunk number 15: iinf.tex:1847-1851
###################################################
hhpdf("iinf-f-dist.pdf", height=7.2, width=9)
tmp <- lapply(qqnorm.dat, histogram, nint=15)
update(do.call(latticeExtra:::c.trellis, c(tmp[-3], list(layout=c(2,3)))), xlab="Observed Values")
hhdev.off()


###################################################
### code chunk number 16: iinf.tex:1869-1890
###################################################
hhpdf("iinf-f1.pdf", height=4, width=7)
points <- ppoints(100)
qqt.data <- data.frame(
  ppoints=points,
  q3=qt(points, 3),
  q5=qt(points, 5),
  q7=qt(points, 7),
  qn=qnorm(points))
qqt.data <- within(qqt.data, {
  d3 <- dt(q3, 3)
  d5 <- dt(q5, 5)
  d7 <- dt(q7, 7)
  dn <- dnorm(qn)
})
head(qqt.data)
xyplot(q3+q5+q7+qn ~ qn, data=qqt.data, type="l", aspect="iso", ylab="distribution",
       par.settings=list(superpose.line=list(lty=c(5,4,2,1))),
       auto.key=list(space="right", border=TRUE,
                     text=c("t, 3 df", "t, 5 df", "t, 7 df", "normal"),
                     lines=TRUE, points=FALSE))
hhdev.off()


###################################################
### code chunk number 17: iinf.tex:1911-1923
###################################################
hhpdf("iinf-f2.pdf", height=4, width=7)
qqt.melt <- cbind(reshape2::melt(qqt.data[,2:5], value.name="quantile",
                                 variable.name="distribution"),
                  reshape2::melt(qqt.data[,9:6], value.name="density"))
levels(qqt.melt$distribution) <- c("t, 3 df", "t, 5 df", "t, 7 df", "normal")
head(qqt.melt)

xyplot(density ~ quantile, group=distribution, data=qqt.melt, type="l",
       par.settings=list(superpose.line=list(lty=c(5,4,2,1))),
       auto.key=list(space="right", border=TRUE, lines=TRUE, points=FALSE)) +
          layer(panel.abline(a=0, b=0, col="gray60"))
hhdev.off()


###################################################
### code chunk number 18: iinf.tex:1986-1993
###################################################
hhcapture("iinf-f3.Rout", '
rt5 <- rt(300, df=5)
rnn <- rnorm(300)

ks.test(rt5, function(x)pt(x, df=2))
ks.test(rnn, function(x)pt(x, df=2))
')


###################################################
### code chunk number 19: iinf.tex:2007-2034
###################################################
plotks1 <- function(x, pfunction, ..., xlim=range(x), ylim.dev=c(-.10, .10)) {
  sx <- sort(x)
  xset <- seq(xlim[1], xlim[2], length=101)
  xref <- xyplot(pfunction(xset) ~ xset, type="l", col="black", xlim=xlim, ylim=c(-.05, 1.05))
  empcdf <- (1:length(sx))/length(sx)
  resizePanels(h=c(1, 2.5),
               update(strip=FALSE, strip.left=TRUE, ylab=NULL, scales=list(rot=0),
                      xlim=xlim, ylim=list(ylim.dev, c(-.05, 1.05)),
               c(deviations=
                   xyplot(empcdf-pfunction(sx) ~ sx, type="h",
                          col=trellis.par.get("superpose.line")$col[1], ...),
                 "cumulative distribution with deviations"=
                   xref + layer(panel.segments(sx, pfunction(sx), sx, empcdf,
                                               col=trellis.par.get("superpose.line")$col[1]),
                           data=list(sx=sx, pfunction=pfunction, empcdf=empcdf)),
                 layout=c(1,2), x.same=TRUE, y.same=FALSE)))
}

KSt5 <- plotks1(rt5, function(x) pt(x, df=2), xlim=range(rt5, rnn),
                xlab="random selection from t with 5 df", main="compare to t with 2 df")
KSnn <- plotks1(rnn, function(x) pt(x, df=2), xlim=range(rt5, rnn),
                xlab="random selection from normal", main="compare to t with 2 df")

hhpdf("iinf-f3.pdf", height=6, width=7)
print(KSt5, position=c(0,0,.5,1), more=TRUE)
print(KSnn, position=c(.5,0,1,1), more=FALSE)
hhdev.off()


###################################################
### code chunk number 20: iinf.tex:2062-2065
###################################################
hhcapture("iinf-f32.Rout", '
ks.test(rt5, rnn)
')


###################################################
### code chunk number 21: iinf.tex:2078-2094
###################################################
plotks2 <- function(x, y, ..., xlim=range(x,y), ylim.dev=c(-.10, .10)) {
  sx <- sort(x)
  sy <- sort(y)
  xset <- seq(xlim[1], xlim[2], length=101)
  xref <- xyplot(pfunction(xset) ~ xset, type="l", col="black", xlim=xlim, ylim=c(-.05, 1.05))
  empcdfx <- (1:length(sx))/length(sx)
  empcdfy <- (1:length(sy))/length(sy)

  xyplot(empcdfx ~ sx, xlim=xlim, ylim=c(-.05, 1.05), type="l", col="red", ...) +
     xyplot(empcdfy ~ sy, xlim=xlim, ylim=c(-.05, 1.05), type="l", col="blue", ...)
}

hhpdf("iinf-f32.pdf", height=3, width=7)
plotks2(rt5, rnn, xlab="random selections from t with 5 df (red) and normal (blue)",
        ylab="empirical CDF", main="Two-Sample Kolmogorov-Smirnov Test")
hhdev.off()


