First example in guided notes

set.seed(12345)

lambda<-2
numsim<-100
mean5<-rep(0,numsim)
mean25<-mean5
mean100<-mean5

for (i in 1:numsim){
  mean5[i]<-mean(rexp(5,lambda))
  mean25[i]<-mean(rexp(25,lambda))
  mean100[i]<-mean(rexp(100,lambda))
  
}
boxplot(mean5,mean25,mean100,names=c("n=5","n=25","n=100"))
title("Distribution of Mean of Exponentials") 

y <- seq(0,1,.1)
#lines(y,(1/sqrt(2*pi))*exp(-.5*(y)^2))

mean(mean5)
## [1] 0.5352591
mean(mean25)
## [1] 0.5027771
mean(mean100)
## [1] 0.4996425
n <- c(5,25,100)
mean <- c(mean(mean5),mean(mean25),mean(mean100))
sd <- c(sd(mean5),sd(mean25),sd(mean100))
kable(cbind(n,mean,sd)) %>%
  kable_paper("hover", full_width = F)
n mean sd
5 0.5352591 0.2551066
25 0.5027771 0.1011640
100 0.4996425 0.0458466

Page 5 Ex

library(kableExtra)
set.seed(12345678)

#Generate from normal(0,1) 3 samples of size 5

s <- 3
n <- 5
x4 <- rep(0,s)
x5 <- rep(0,s)
x6 <- rep(0,s)

mymatrix <- matrix(rnorm(s*n), nrow=s,ncol=n)
mydata <- cbind(mymatrix,x4,x5,x6)

mydata[,6] <- apply(mydata[,1:5],1,mean)
mydata[,7] <- apply(mydata[,1:5],1,function(s){mean(s,trim=.2)} )
mydata[,8] <- apply(mydata[,1:5],1,median)

thetaMC <- rep(0,3)
biasMC<- rep(0,3)
SDMC<- rep(0,3)
MSEMC<- rep(0,3)



for (i in 1:3) {
  thetaMC[i] <- mean(mydata[,n+i])
  biasMC[i] <- thetaMC[i] - 0
  SDMC[i] <- sd(mydata[,n+i])
  MSEMC[i] <- (SDMC[i])^2 + (biasMC[i])^2
}
estimator <- c("sample mean", "trimmed mean", "median")

kable(round(mydata,3)) %>%
  kable_paper("hover", full_width = F)
x4 x5 x6
1.148 -0.257 -0.756 2.552 -0.280 0.481 0.204 -0.257
1.618 0.892 0.984 -0.525 -0.052 0.583 0.608 0.892
-2.379 -0.752 1.099 -0.221 -0.497 -0.550 -0.490 -0.497 
kable(cbind(estimator,round(thetaMC,3),round(biasMC,3),round(SDMC,3),round(MSEMC,3)))%>%
  kable_styling()
estimator
sample mean 0.172 0.172 0.627 0.423
trimmed mean 0.107 0.107 0.555 0.32
median 0.046 0.046 0.742 0.553 
mu <- 0
m <- 1000
g <- numeric(m)
MSE.Num<-numeric(m)
    
for (i in 1:m) {
  x <- rnorm(2)
  g[i] <- abs(x[1] - x[2])
  MSE.Num[i]<-(g[i]-mu)^2 # MSE calculation of the numerator
}
    est <- mean(g)
    est
## [1] 1.094918
    bias.g<-est-0 # Monte Carlo Bias, note mu ==0

 bias.g
## [1] 1.094918
MC.SE<-sd(g) # MC Standard error
var.theor<-2-(4/pi) # theoretical variance
MC.SE
## [1] 0.8096886
var.theor
## [1] 0.7267605
MC.MSE<-sum(MSE.Num)/m # Monte Carlo MSE
MC.MSE.shortcut<-(bias.g)^2+(MC.SE)^2 # Short cut way to calculate MSE

MC.MSE
## [1] 1.853786
MC.MSE.shortcut
## [1] 1.854442
names <- c("mean estimate","bias", "SE","theor. variance", "MSE")
kable(rbind(c("mean estimate","bias", "SE","theor. variance", "MSE"), round(c(est,bias.g,MC.SE,var.theor,MC.MSE),3))) %>%
  kable_styling()
mean estimate bias SE theor. variance MSE
1.095 1.095 0.81 0.727 1.854

6.4 and 6.5

 n <- 20
    alpha <- .05
    UCL <- replicate(1000, expr = {
        x <- rnorm(n, mean = 0, sd = 2)
        (n-1) * var(x) / qchisq(alpha, df = n-1)
        } )
    #count the number of intervals that contain sigma^2=4
    sum(UCL > 4)
## [1] 951
    #or compute the mean to get the confidence level
    mean(UCL > 4)
## [1] 0.951

Page 11 ex

set.seed(1234)
n <- 10
alpha <- .025
UCL <- replicate(10000, expr = {
   x <- rnorm(n, mean = 2, sd = 1) +
   (sd(x)/sqrt(n)) * qt(p=(1-alpha/2), df = (n-1))
    } )

sum(UCL > 2)    
## [1] 52584
mean(UCL > 2)
## [1] 0.52584

Ex 6.7 pg 163 in text

set.seed(12345)
n <- 20
alpha <- .05
mu0 <- 500
sigma <- 100

m <- 10000          #number of replicates
p <- numeric(m)     #storage for p-values
  for (j in 1:m) {
     x <- rnorm(n, mu0, sigma)
     ttest <- t.test(x, alternative = "greater", mu = mu0)
        p[j] <- ttest$p.value
        }

p.hat <- mean(p < alpha)
se.hat <- sqrt(p.hat * (1 - p.hat) / m)
kable(rbind(c("p.hat","se.hat"),c(round(p.hat,3), round(se.hat,3)))) %>%
  kable_paper(full_width = F)
p.hat se.hat
0.051 0.002

Different version of 6.7 (Chi squared)

set.seed(23456)
n2 <- 20
alpha2 <- .05
mu02 <- 500


m2 <- 10000          #number of replicates
q <- numeric(m2)     #storage for p-values
  for (j in 1:m2) {
     x2 <- rchisq(n2, df = 500)
     ttest2 <- t.test(x2, alternative = "greater", mu = mu02)
        q[j] <- ttest2$p.value
        }

p.hat2 <- mean(q < alpha)
se.hat2 <- sqrt(p.hat2 * (1 - p.hat2) / m2)
kable(rbind(c("p.hat2","se.hat2"),c(round(p.hat2,3), round(se.hat2,3)))) %>%
  kable_paper(full_width = F)
p.hat2 se.hat2
0.044 0.002

Power curve for ex 6.7

library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
set.seed(34567)
n3 <- 20
mu03 <- 500
s <- 10000
sigma3 <- 100
mu <- c(seq(450, 650, 10))
M <- length(mu)
power <- numeric(M)

r <- numeric(s)     #storage for p-values
  for (i in 1:M) {
        mu1 <- mu[i]
        pvalues <- replicate(s, expr = {
            #simulate under alternative mu1
            x <- rnorm(n3, mean = mu1, sd = sigma3)
            ttest3 <- t.test(x, alternative = "greater", mu = mu03)
            ttest3$p.value } )
        power[i] <- mean(pvalues <= .05)
  }
plot(mu,power)
abline(v = mu03, lty = 1)
abline(h = .05, lty = 1)

se <- sqrt(power * (1-power) / m)
    errbar(mu, power, yplus = power+se, yminus = power-se,
        xlab = bquote(theta))
    lines(mu, power, lty=3)

    detach(package:Hmisc)
    par(ask = FALSE)