n.ajstu <- 100000
mean.ajstu <- 100
sd.ajstu <- 10

set.seed(1024)
ajstu <- rnorm2(n.ajstu, mean=mean.ajstu, sd=sd.ajstu)

mean(ajstu)
sd(ajstu)
var(ajstu)

iter <- 10000 # # of sampling 

n.4 <- 4
means4 <- rep (NA, iter)
for(i in 1:iter){
  means4[i] = mean(sample(ajstu, n.4))
}

n.25 <- 25
means25 <- rep (NA, iter)
for(i in 1:iter){
  means25[i] = mean(sample(ajstu, n.25))
}

n.100 <- 100
means100 <- rep (NA, iter)
for(i in 1:iter){
  means100[i] = mean(sample(ajstu, n.100))
}

n.400 <- 400
means400 <- rep (NA, iter)
for(i in 1:iter){
  means400[i] = mean(sample(ajstu, n.400))
}

n.900 <- 900
means900 <- rep (NA, iter)
for(i in 1:iter){
  means900[i] = mean(sample(ajstu, n.900))
}

n.1600 <- 1600
means1600 <- rep (NA, iter)
for(i in 1:iter){
  means1600[i] = mean(sample(ajstu, n.1600))
}

n.2500 <- 2500
means2500 <- rep (NA, iter)
for(i in 1:iter){
  means2500[i] = mean(sample(ajstu, n.2500))
}

h4 <- hist(means4)
h25 <- hist(means25)
h100 <- hist(means100)
h400 <- hist(means400)
h900 <- hist(means900)
h1600 <- hist(means1600)
h2500 <- hist(means2500)


plot(h4, ylim=c(0,3000), col="red")
plot(h25, add = T, col="blue")
plot(h100, add = T, col="green")
plot(h400, add = T, col="grey")
plot(h900, add = T, col="yellow")

m4 <- mean(means4)
m25 <- mean(means25)
m100 <- mean(means100)
m400 <- mean(means400)
m900 <- mean(means900)
m1600 <- mean(means1600)
m2500 <- mean(means2500)

s4 <- sd(means4)
s25 <- sd(means25)
s100 <- sd(means100)
s400 <- sd(means400)
s900 <- sd(means900)
s1600 <- sd(means1600)
s2500 <- sd(means2500)

sss <- c(4,25,100,400,900,1600,2500) # sss sample sizes
means <- c(m4, m25, m100, m400, m900, m1600, m2500)
sds <- c(s4, s25, s100, s400, s900, s1600, s2500)

temp <- data.frame(sss,
                   means,
                   sds)

temp

ses <- rep (NA, length(sss)) # std error memory
for(i in 1:length(sss)){
  ses[i] = sqrt(var(ajstu)/sss[i])  # std errors by theorem
}

data.frame(ses)
se.1 <- ses
se.2 <- 2 * ses 

lower.s2 <- mean(ajstu)-se.2
upper.s2 <- mean(ajstu)+se.2
data.frame(cbind(sss, ses, lower.s2, upper.s2))

# 12/2 lecture 
# note that we draw the statistical calculation 
# by "diff/se" = "diff/random_error"
n <- 80
mean.sample <- 103

sample <- rnorm2(n, mean.sample, sd.ajstu)
mean(sample)
sd(sample)

diff <- mean.sample - mean.ajstu # this is actual difference
se <- sd.ajstu / sqrt(n) # this is random error 
t.cal <- diff/se
t.cal
qnorm(0.025, lower.tail = F)
qnorm(0.01/2, lower.tail = F)
qt(0.05/2, n-1, lower.tail=F)

t.test(sample, mu=mean.ajstu)

# or we obtain the exact p value
p.value <- pt(t.cal, n-1, lower.tail = F)
p.value*2