COMMunication
RESearch.NET

This is an old revision of the document!

---

# variance 
# why we use n-1 when calculating 
# population variance (instead of N)
a <- rnorm2(100000000, 100, 10)
a.mean <- mean(a)
ss <- sum((a-a.mean)^2)
n <- length(a)
df <- n-1
ss/n
ss/df

# standard deviation 68, 85, 99 rule 
one.sd <- .68
two.sd <- .95
thr.sd <- .99
1-one.sd
(1-one.sd)/2
qnorm(0.16)
qnorm(1-0.16)

1-two.sd
(1-two.sd)/2
qnorm(0.025)
qnorm(0.975)

1-thr.sd
(1-thr.sd)/2
qnorm(0.005)
qnorm(1-0.005)


# sampling distribution 
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)
min.value <- min(ajstu)
max.value <- max(ajstu)
min.value
max.value

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")


sss <- c(4,25,100,400,900,1600,2500) # sss sample sizes
ses <- rep (NA, length(sss)) # std errors
for(i in 1:length(sss)){
  ses[i] = sqrt(var(ajstu)/sss[i])
}
ses.means4 <- sqrt(var(means4))
ses.means25 <- sqrt(var(means25))
ses.means100 <- sqrt(var(means100))
ses.means400 <- sqrt(var(means400))
ses.means900 <- sqrt(var(means900))
ses.means1600 <- sqrt(var(means1600))
ses.means2500 <- sqrt(var(means2500))
ses.real <- c(ses.means4, ses.means25,
              ses.means100, ses.means400,
              ses.means900, ses.means1600,
              ses.means2500)
ses.real 

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, ses.real, lower.s2, upper.s2))

min.value <- min(ajstu)
max.value <- max(ajstu)
min.value
max.value
# means25 분포에서 population의 가장 최소값인 
# min.value가 나올 확률은 얼마나 될까?
m.25 <- mean(means25)
sd.25 <- sd(means25)
m.25
sd.25
pnorm(min.value, m.25, sd.25)
# 위처럼 최소값, 최대값의 probability를 가지고
# 확률을 구하는 것은 의미가 없음. 즉, 어떤 샘플의
# 평균이 원래 population에서 나왔는지를 따지는 것을
# 그 pop의 최소값으로 판단하는 것은 의미가 없음.
# 그 보다는 sd.25와 같은 standard deviation을 가지
# 고 보는 것이 타당함. 이 standard deviation 을
# standard error라고 부름
se.25 <- sd.25
se.25.2 <- 2 * sd.25
mean(means25)+c(-se.25.2, se.25.2)
# 위처럼 95% certainty를 가지고 구간을 판단하거나
# 아래 처럼 점수를 가지고 probability를 볼 수도
# 있음
pnorm(95.5, m.25, sd.25)
# 혹은 
pnorm(95.5, m.25, sd.25) * 2
# 매우중요. 위의 이야기는 내가 25명의 샘플을 취
# 했을 때, 그 점수가 95.5가 나올 확률을 구해 본것
# 다른 예로
pnorm(106, m.25, sd.25, lower.tail = F)
# 혹은 
2 * pnorm(106, m.25, sd.25, lower.tail = F)