## Chapter 6
## November 17, 2007
## Alan T. Arnholt
#
options(width=65)
library(PASWR)
# Example 6.1
Combinations(5,3)
#
# Example 6.2
t(SRS(c(2,5,8,12,13), 3))
#
# Example 6.3
sample(1:180, 5, replace=FALSE)
#
# Example 6.4
sample(x=(1:20),size=5,prob=c(rep(1/26,18),rep(4/26,2)),replace=FALSE)
#
# 1 in 100 Systematic Sampling strategy
seq(sample(1:100,1), 1000, 100)
#
# Example 6.6
seq(sample(1:50,1), 1000, 50)
#
# Example 6.8
rolls <- sample(1:6,100, replace=TRUE)
table(rolls)
table(rolls)/100 # epdf
# similar to Figure 6.1
plot(ecdf(rolls))
#
# Example 6.10
N <- 6
n <- 2
pop <- 1:N
rs <- expand.grid(Draw1=pop, Draw2=pop) # Possible random samples
xbarN <- apply(rs, 1, mean) # Means of all rs values
s2N <- apply(rs, 1, var)    # Variance of all rs values
TOT1<- cbind(rs, xbarN=xbarN, s2N=s2N)
TOT1            # Numerical values for Table 6.5
table(xbarN)    # Numerical values for Table 6.6
table(s2N)      # Numerical values for Table 6.7
MU <- mean(pop) # Population mean
VAR <- sum((pop-mean(pop))^2)*(1/N) # Population variance
MU.xbarN <- mean(xbarN)             # Expected value of xbarN
E.s2N <- mean(s2N)                  # Expected value of s2N
VAR.xbarN <- sum((xbarN-mean(xbarN))^2)*(1/(N*N))
reN <- c(MU, MU.xbarN, VAR, E.s2N, VAR.xbarN)
names(reN)<-c("MU", "MU.xbarN", "VAR", "E.s2N", "V.xbarN")
reN                 # Numerical values for Case 1 in Table 6.12
srs <- SRS(1:N, n)  # Possible simple random samples
xbari <- apply(srs, 1, mean)  # Means of simple random samples
s2i <- apply(srs, 1, var)     # Variances of simple random samples
TOT <- cbind(srs, xbari, s2i)
dimnames(TOT)[[2]]<-c("Draw1", "Draw2", "xbari", "s2i")
TOT               # Numerical values for Table 6.8
table(xbari)      # Numerical values for Table 6.9
table(s2i)        # Numerical values for Table 6.10
MU <- mean(pop)   # Population mean
VAR <- sum((pop-mean(pop))^2)*(1/N)   # Population variance
MU.xbar <- mean(xbari)                # Expected value of xbari
E.s2 <- mean(s2i)                     # Expected value of s2i
VAR.xbar <- sum((xbari-mean(xbari))^2)*(1/choose(N,n))
results <- c(MU, MU.xbar, VAR, E.s2, VAR.xbar)
names(results)<-c("MU", "MU.xbari", "VAR", "E.s2i", "V.xbari")
print(results)    # Numerical values for Case 2 in Table 6.12
#
#
# Example 6.11
nsize(b=3,sigma=12)
#
# Example 6.12
pnorm(50,50,1.897367)
pnorm(75,50,1.897367)
pnorm(100,50,1.897367)
#
# Example 6.13
# part c.
nsize(b=1,sigma=2)
#
#
################################################################################
#  n2UNIFsim code                                                              #
################################################################################


par(mfrow=c(2,2))
set.seed(15)
m <- 50000
n <-2
lv <- (2-(2+sqrt(12))/2)/(1/5)
uv <- (2+sqrt(12))/(2/5)
mu <- (lv+uv)/2
SIGMA <- (uv-lv)/sqrt(12)
MEANSunif <- array(0,m)
for (i in 1:m) {MEANSunif[i]<-mean(runif(n,lv,uv))}
MEANSexp <- array(0,m)
for (i in 1:m) {MEANSexp[i]<-mean(rexp(n,1/5))}
ll <- mu - 5*SIGMA/sqrt(n)
ul <- mu + 5*SIGMA/sqrt(n)
hist(MEANSunif,prob=T,xlab="xbar",nclass="scott", col=2,xlim=c(ll,ul),
ylim=c(0,1.5*dnorm(mu,mu,SIGMA/sqrt(n))),main="A")
lines(seq(ll,ul,.001), dnorm(seq(ll,ul,.001),mu,SIGMA/sqrt(n)),lwd=2)
Mun <- mean(MEANSunif)
sun <- sqrt(var(MEANSunif))
UNIFpercents <- c(sum(MEANSunif>Mun-50*sun & MEANSunif<Mun-2*sun),
sum(MEANSunif>Mun-2*sun & MEANSunif<Mun-1*sun),
sum(MEANSunif>Mun-1*sun & MEANSunif<Mun),
sum(MEANSunif>Mun & MEANSunif<Mun+1*sun),
sum(MEANSunif>Mun+1*sun & MEANSunif<Mun+2*sun),
sum(MEANSunif>Mun+2*sun & MEANSunif<Mun+50*sun) ) / m
hist(MEANSexp,prob=T,xlab="xbar",nclass="scott", col=2,xlim=c(ll,ul),
ylim=c(0,1.5*dnorm(mu,mu,SIGMA/sqrt(n))),main="B")
lines(seq(ll,ul,.001), dnorm(seq(ll,ul,.001),mu,mu/sqrt(n)),lwd=2)
Mex <- mean(MEANSexp)
sex <- sqrt(var(MEANSexp))
EXPpercents <- c( sum(MEANSexp>Mex-50*sex & MEANSexp<Mex-2*sex),
sum(MEANSexp>Mex-2*sex & MEANSexp<Mex-1*sex),
sum(MEANSexp>Mex-1*sex & MEANSexp<Mex),
sum(MEANSexp>Mex & MEANSexp<Mex+1*sex),
sum(MEANSexp>Mex+1*sex & MEANSexp<Mex+2*sex),
sum(MEANSexp>Mex+2*sex & MEANSexp<Mex+50*sex) ) / m
print(c(Mun,sun))
print(c(Mex,sex))
print(UNIFpercents)
print(EXPpercents)
#
#
set.seed(15)
m <- 50000
n <-16
lv <- (2-(2+sqrt(12))/2)/(1/5)
uv <- (2+sqrt(12))/(2/5)
mu <- (lv+uv)/2
SIGMA <- (uv-lv)/sqrt(12)
MEANSunif <- array(0,m)
for (i in 1:m) {MEANSunif[i]<-mean(runif(n,lv,uv))}
MEANSexp <- array(0,m)
for (i in 1:m) {MEANSexp[i]<-mean(rexp(n,1/5))}
ll <- mu - 5*SIGMA/sqrt(n)
ul <- mu + 5*SIGMA/sqrt(n)
hist(MEANSunif,prob=T,xlab="xbar",nclass="scott", col=2,xlim=c(ll,ul),
ylim=c(0,1.5*dnorm(mu,mu,SIGMA/sqrt(n))),main="C")
lines(seq(ll,ul,.001), dnorm(seq(ll,ul,.001),mu,SIGMA/sqrt(n)),lwd=2)
Mun <- mean(MEANSunif)
sun <- sqrt(var(MEANSunif))
UNIFpercents <- c(sum(MEANSunif>Mun-50*sun & MEANSunif<Mun-2*sun),
sum(MEANSunif>Mun-2*sun & MEANSunif<Mun-1*sun),
sum(MEANSunif>Mun-1*sun & MEANSunif<Mun),
sum(MEANSunif>Mun & MEANSunif<Mun+1*sun),
sum(MEANSunif>Mun+1*sun & MEANSunif<Mun+2*sun),
sum(MEANSunif>Mun+2*sun & MEANSunif<Mun+50*sun) ) / m
hist(MEANSexp,prob=T,xlab="xbar",nclass="scott", col=2,xlim=c(ll,ul),ylim=c(0,1.5*dnorm(mu,mu,SIGMA/sqrt(n))),main="D")
lines(seq(ll,ul,.001), dnorm(seq(ll,ul,.001),mu,mu/sqrt(n)),lwd=2)
Mex <- mean(MEANSexp)
sex <- sqrt(var(MEANSexp))
EXPpercents <- c( sum(MEANSexp>Mex-50*sex & MEANSexp<Mex-2*sex),
sum(MEANSexp>Mex-2*sex & MEANSexp<Mex-1*sex),
sum(MEANSexp>Mex-1*sex & MEANSexp<Mex),
sum(MEANSexp>Mex & MEANSexp<Mex+1*sex),
sum(MEANSexp>Mex+1*sex & MEANSexp<Mex+2*sex),
sum(MEANSexp>Mex+2*sex & MEANSexp<Mex+50*sex) ) / m
print(c(Mun,sun))
print(c(Mex,sex))
print(UNIFpercents)
print(EXPpercents)
par(mfrow=c(1,1))
# Name running n=2 and n=16 n2unifexp4.ps
# Name running n=36 and n=100 n36unifexp4.ps
#
################################################################################
#
# Example 6.14
round(1 - pnorm(4.5,4,sqrt(0.03)), 3)
#
# Example 6.16
set.seed(4)
m <- 1000
nx <- 100
ny <- 81
mux <- 100
sigx <- 10
muy <- 50
sigy <- 9
muxy <- mux - muy
sigxy <- sqrt((sigx^2/nx) + (sigy^2/ny))
meansX <- array(0, m) # Array of m zeros
meansY <- array(0, m) # Array of m zeros
for(i in 1:m) {meansX[i] <- mean(rnorm(nx, mux, sigx))}
for(i in 1:m) {meansY[i] <- mean(rnorm(ny, muy, sigy))}
XY <- meansX - meansY
ll <- muxy - 3.4 * sigxy
ul <- muxy + 3.4 * sigxy
hist(XY, prob = TRUE, xlab = "xbar-ybar", nclass = "scott",col = 13,
xlim = c(ll, ul), ylim = c(0, 0.3),main="",ylab="")
# name: xbar.ybar.ps
lines(seq(ll, ul, 0.05),dnorm(seq(ll, ul, 0.05),muxy,sigxy),lwd = 3)
print(round(c(mean(XY), sqrt(var(XY))), 2))
sum(XY < 52)/1000
round(pnorm(52, 50, sqrt(2)), 2)
#
# Example 6.17
# part a.
1 - pnorm(0.11,0.10,sqrt(0.1*0.9/500))
1 - pnorm(55,500*0.1,sqrt(500*0.1*0.9))
1 - pbinom(54,500,0.10)
# part b.
pnorm(0.12,0.10,sqrt(0.1*0.9/500))
pnorm(60,500*.10,sqrt(500*0.1*0.9))
pbinom(60,500,0.1)
#
# Example 6.18
# Approximate answer
sig <- sqrt((0.383*0.617)/250)
pnorm(0.402,0.383,sig) - pnorm(0.358,0.383,sig)
# Exact answer
pbinom(100,250,.383) - pbinom(89,250,.383)
# Or
sum(dbinom(90:100,250,.383))
#
# Example 6.19
# part a.
1 - pchisq(126,150)
# part b.
pchisq(50,65) - pchisq(40,65)
# part c.
1 - pchisq(260,220)
# part d.
qchisq(.6,100)
#
# Example 6.20
# part a.
pchisq(24,10)
# part b.
1 - pchisq(4.87,10)
# part c.
qchisq(0.50,10)
#
# Example 6.22
pchisq(15.99,10) - pchisq(4.87,10)
#
# Example 6.23
1 - pchisq(43.78,30)
#
################################################################################
# Example 6.24
set.seed(302)
par(mfrow=c(2,2))
m <- 1000; n <- 15
varNC14 <- array(0,m)         # Array with m zeros
for (i in 1:m) {varNC14[i] <- var(rnorm(n))}
NC14 <- (n-1)*varNC14/1
hist(NC14,prob=TRUE,ylim=c(0,0.09),xlab="NC14",col=2,xlim=c(0,60),
nclass="scott", main="", ylab="")
lines(seq(0,60,.1), dchisq(seq(0,60,.1), n-1), lwd=3)
varEC14 <- array(0,m)
for (i in 1:m) {varEC14[i] <- var(rexp(n))}
EC14 <- (n-1)*varEC14/1
hist(EC14,prob=TRUE,ylim=c(0,0.09),xlab="EC14",col=4,xlim=c(0,60),
nclass="scott", main="", ylab="")
lines(seq(0,60,.1), dchisq(seq(0,60,.1), n-1), lwd=3)
n <- 100
varNC99 <- array(0,m)
for (i in 1:m) {varNC99[i] <- var(rnorm(n))}
NC99 <- (n-1)*varNC99/1
hist(NC99,prob=TRUE,ylim=c(0,0.03),xlab="NC99",col=2,xlim=c(0,210),
nclass="scott", main="", ylab="")
lines(seq(0,210,.1), dchisq(seq(0,210,.1), n-1), lwd=3)
varEC99 <- array(0,m)
for (i in 1:m) {varEC99[i] <- var(rexp(n))}
EC99 <- (n-1)*varEC99/1
hist(EC99,prob=TRUE,ylim=c(0,0.03),xlab="EC99",col=4,xlim=c(0,210),
nclass="scott", main="", ylab="")
lines(seq(0,210,.1), dchisq(seq(0,210,.1), n-1), lwd=3)
NC14 <- c(mean(varNC14), var(varNC14), mean(NC14), var(NC14))
EC14 <- c(mean(varEC14), var(varEC14), mean(EC14), var(EC14))
NC99 <- c(mean(varNC99), var(varNC99), mean(NC99), var(NC99))
EC99 <-c (mean(varEC99), var(varEC99), mean(EC99), var(EC99))
MAT <- round(rbind(NC14, EC14, NC99, EC99), 4)
colNAM <- c("E(S^2)", "Var(S^2)", "E(X^2)", "Var(X^2)")
rowNAM <- c("NC14", "EC14" ,"NC99", "EC99")
dimnames(MAT) <- list(rowNAM ,colNAM)
print(MAT)          # Numerical values for Table 6.14
# Name: SIMCHI2.ps
par(mfrow=c(1,1))
#
# Example 6.26
qf(0.95,5,10)
qf(0.05,5,10)
#
# Example 6.27
qf(.975,19,19)
qf(.025,19,19)
################################################################################
################################################################################