COMMunication
RESearch.NET

This is an old revision of the document!

---

n.ajstu <- 100000
mean.ajstu <- 110
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))
}

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

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

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")
# plot(h2500, add = T, col="brown")
# plot(h3600, add = T, col="orange")

sss <- c(4,25,100,400,900,1600,2500,3600) # sss sample sizes
ses <- rep (NA, length(sss)) # std errors
for(i in 1:length(sss)){
  ses[i] = sqrt(var(ajstu)/sss[i])
}

data.frame(ses)
se.1 <- ses
se.2 <- 2 * ses 
# alt.2 <- qnorm(.975)
# se.2 <- alt.2 * ses 
se.2

lower.s2 <- mean(ajstu)-se.2
upper.s2 <- mean(ajstu)+se.2

sample <- rnorm2(100, 112.41, 10)


ses <- round(ses,3)
smean <- mean(sample)
diff <- smean - mean(ajstu)
diff 

# z-test 
z.sc <- (diff)/ses
z.sc
# z.crit <- 2
z.crit <- qnorm(.975)
z.crit
data.frame(cbind(sss, smean, mean(ajstu), diff, ses, 
                 lower.s2, upper.s2, 
                 z.sc, z.crit, z.sc>z.crit))

# t-test
ses <- round(ses,3)
smean <- smean
diff <- smean - mean(ajstu)
dfs <- sss - 1
t.sc <- diff/ses
perc <- round(pt(t.sc, dfs, lower.tail = FALSE),9)
t.crit <- qt(.975, dfs, lower.tail = TRUE)
data.frame(cbind(sss, smean, mean(ajstu), diff, ses, 
                 t.sc, dfs, perc, t.crit, t.sc>t.crit))

# 참고
qt(.975, 10000000000, lower.tail = TRUE)
qnorm(.975)

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

sa <- rnorm2(100, 110, 10)
sb <- rnorm2(100, 115, 10)

sa.mean <- mean(sa)
sb.mean <- mean(sb)

na <- length(sa)
nb <- length(sb)
dfa <- na-1
dfb <- nb-1

var(sa)
var(sb)

ssa <- var(sa)*dfa
ssb <- var(sb)*dfb
ssa
ssb

ssa+ssb
dfa+dfb

var.pooled <- (ssa+ssb)/(dfa+dfb)
vp <- var.pooled
vp

se <- sqrt((vp/na)+(vp/nb))
diff <- sa.mean - sb.mean 
t.cal <- diff/se
se
diff
t.cal
t.crit <- qt(.975, 198)
t.crit
pt(t.cal, 198)*2

t.test(sa, sb)