COMMunication
RESearch.NET

n.ajstu <- 100000 # 모집단 100000
mean.ajstu <- 110 # 모집단 평균 110이라 가정
sd.ajstu <- 10 # 표준편차 10이라 가정

set.seed(1024)
# rnorm2 펑션을 이용해서 모집단 만들기 
rnorm2 <- function(n,mean,sd) { mean+sd*scale(rnorm(n)) } 
ajstu <- rnorm2(n.ajstu, mean=mean.ajstu, sd=sd.ajstu)

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


iter <- 10000 # # of sampling 

# 이 후 n.# 은 샘플사이즈를 말함
n.4 <- 4   # 샘플사이즈 4
means4 <- rep (NA, iter)
# sample 평션을 이용해서 4개 샘플을 ajoust에서
# 취한 후, 평균을 내서 이를 means[i]에 기록한다
# 이것을 iteration 숫자만큼 한다
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")
# plot(h2500, add = T, col="brown")

# 아래는 종합적인 출력을 위해서 좀 복잡하게 
# 코딩을 하였지만, 요지는 특정 샘플사이즈에서 
# standard deviation of sample means 즉 
# standard error를 구하는 과정
sss <- c(4,25,100,400,900,1600,2500,3600) # 8 종류의 샘플사이즈
ses <- rep (NA, length(sss)) # 각 샘플 사이즈에 대응하는 standard error 값 저장을 위한 공간
for(i in 1:length(sss)){
  ses[i] = sqrt(var(ajstu)/sss[i])
}

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

# 2 standard deviation of sample means를 (se) 썼을 때 
# 아랫쪽, 위쪽 경계
lower.s2 <- mean(ajstu)-se.2
upper.s2 <- mean(ajstu)+se.2


# 이제 머리가 좋아지는 약을 먹은 하나의 샘플을
# 가정하고 이 샘플의 평균이  113라고 가정하면 
smean <- 113

# z-test를 위한 계산에서 분자는 difference라고 했고 
# 이를 구한 것
diff <- smean - mean(ajstu)
diff 

# z-test diff를 stanard error 값으로 나눈 것 
# 이 diff를 샘플사이즈 별로 (8종류) 구한 ses로 
# 나누어 본다
z.sc <- (diff)/ses
z.sc
# 위의 z.sc 점수와 강사가 이야기한 2점과 비교
# 이 2점을 (standard error 2개) z critical score라고
# 부른다
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 의 상황에서는 z critical value가 2였던 것에
# 비해서 샘플사이즈에 따라서 critical value가 달라지게
# 된다. 이를 구한 것이 아래의 qt평션을 이용한 t.crit 
diff <- smean - mean(ajstu)
dfs <- sss - 1
t.sc <- diff/ses
# 아래는 샘플사이즈 별로 각각 구한 
# t calculated value 들에 해당하는 
# t distribution에서의 퍼센트  
perc <- round(pt(t.sc, dfs, lower.tail = FALSE), 5)

# 아래는 degrees of freedom이 dfs 이고 
# 97.5 퍼센트일 때 봐야 할 (구해야 할) 
# t 값을 구한 것. 이것을 샘플사이즈 별로 한 것
# 4, 25, . . . .
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)

# 하나만 따로 생각해 봄
# 평균 113, 표준편차가 10인 100개로 이루어진 
# 샘플을 구하여 sa 에 지정
sa <- rnorm2(100, 113, 10)

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

# two sample t test 의 경우
sa <- rnorm2(100, 110, 10)
sb <- rnorm2(100, 113, 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)