====== Sampling Distribution and z-test ======
<code>
rm(list=ls())
rnorm2 <- function(n,mean,sd){ 
  mean+sd*scale(rnorm(n)) 
}

se <- function(sample) {
  sd(sample)/sqrt(length(sample))
}

ss <- function(x) {
  sum((x-mean(x))^2)
}

################################
N.p <- 1000000
m.p <- 100
sd.p <- 10

p1 <- rnorm2(N.p, m.p, sd.p)
mean(p1)
sd(p1)

p2 <- rnorm2(N.p, m.p+5, sd.p)
mean(p2)
sd(p2)

m.p1 <- mean(p1)
sd.p1 <- sd(p1)
var(p1)


hist(p1, breaks=50, col = rgb(1, 1, 1, 0.5),
     main = "histogram of p1 and p2",)
abline(v=mean(p1), col="black", lwd=3)
hist(p2, add=T, breaks=50, col=rgb(1,1,.5,.5))
abline(v=mean(p2), col="red", lwd=3)


hist(p1, breaks=50, col=rgb(0,.5,.5,.5))
abline(v=mean(p1),lwd=2)
abline(v=mean(p1)-sd(p1), lwd=2)
abline(v=mean(p1)+sd(p1), lwd=2)
abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
abline(v=c(m.p1-3*sd.p1, m.p1+3*sd.p1), lwd=2, col='green')

# area bet black = 68%
# between red = 95%
# between green = 99%

pnorm(m.p1+sd.p1, m.p1, sd.p1)
pnorm(m.p1-sd.p1, m.p1, sd.p1)
pnorm(m.p1+sd.p1, m.p1, sd.p1) - 
  pnorm(m.p1-sd.p1, m.p1, sd.p1)

pnorm(m.p1+2*sd.p1, m.p1, sd.p1) - 
  pnorm(m.p1-2*sd.p1, m.p1, sd.p1)

pnorm(m.p1+3*sd.p1, m.p1, sd.p1) - 
  pnorm(m.p1-3*sd.p1, m.p1, sd.p1)

m.p1
sd.p1

(m.p1-m.p1)/sd.p1
((m.p1-sd.p1) - m.p1) / sd.p1
(120-100)/10 
pnorm(1)-pnorm(-1)
pnorm(2)-pnorm(-2)
pnorm(3)-pnorm(3)

1-pnorm(-2)*2
pnorm(2)-pnorm(-2)

pnorm(120, 100, 10)
pnorm(2)-pnorm(-2)

zscore <- (120-100)/10
pnorm(zscore)-pnorm(-zscore)
zscore

# reminder. 
pnorm(-1)
pnorm(1, 0, 1, lower.tail = F)
pnorm(110, 100, 10, lower.tail = F)
zscore <- (110-100)/10
pnorm(zscore, lower.tail = F)

pnorm(118, 100, 10, lower.tail = F)
pnorm(18/10, lower.tail = F)

z.p1 <- (p1-mean(p1))/sd(p1)
mean(z.p1)
round(mean(z.p1),10)
sd(z.p1)
pnorm(1.8)-pnorm(-1.8)

hist(z.p1, breaks=50, col=rgb(1,0,0,0))
abline(v=c(m.p1, -1.8, 1.8), col='red')
1-(pnorm(1.8)-pnorm(-1.8))


pnorm(1)-pnorm(-1)
1-(pnorm(-1)*2)

pnorm(2)-pnorm(-2)
1-(pnorm(-2)*2)

1-(pnorm(-3)*2)

#
hist(p1, breaks=50, col=rgb(.9,.9,.9,.9))
abline(v=mean(p1),lwd=2)
abline(v=mean(p1)-sd(p1), lwd=2)
abline(v=mean(p1)+sd(p1), lwd=2)
abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
abline(v=c(m.p1-3*sd.p1, m.p1+3*sd.p1), lwd=2, col='green')

# 68%
a <- qnorm(.32/2)
b <- qnorm(1-.32/2)
c(a, b)
c(-1, 1)
# note that
.32/2
pnorm(-1)
qnorm(.32/2)
qnorm(pnorm(-1))

# 95%
c <- qnorm(.05/2)
d <- qnorm(1-.05/2)
c(c, d)
c(-2,2)

# 99% 
e <- qnorm(.01/2)
f <- qnorm(1-.01/2)
c(e,f)
c(-3,3)


pnorm(b)-pnorm(a)
c(a, b)
pnorm(d)-pnorm(c)
c(c, d)
pnorm(f)-pnorm(e)
c(e, f)

qnorm(.5)
qnorm(1)


################################
s.size <- 50

means.temp <- c()
s1 <- sample(p1, s.size, replace = T)
mean(s1)
means.temp <- append(means.temp, mean(s1))
means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
means.temp

iter <- 1000000
# means <- c()
means <- rep(NA, iter)
for (i in 1:iter) {
  # means <- append(means, mean(sample(p1, s.size, replace = T)))
  means[i] <- mean(sample(p1, s.size, replace = T))
}
length(means)
mean(means)
sd(means)
se.s <- sd(means)

hist(means, breaks=50, 
     xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)), 
     col=rgb(1, 1, 1, .5))
abline(v=mean(means), col="black", lwd=3)
# now we want to get sd of this distribution
lo1 <- mean(means)-se.s
hi1 <- mean(means)+se.s
lo2 <- mean(means)-2*se.s
hi2 <- mean(means)+2*se.s
lo3 <- mean(means)-3*se.s
hi3 <- mean(means)+3*se.s
abline(v=mean(means), col="black", lwd=2)
# abline(v=mean(p2), colo='darkgreen', lwd=2)
abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3), 
       col=c("red","red", "blue", "blue", "orange", "orange"), 
       lwd=2)

# meanwhile . . . .
se.s

se.z <- sqrt(var(p1)/s.size)
se.z <- c(se.z)
se.z

se.z.adjusted <- sqrt(var(s1)/s.size)
se.z.adjusted

# sd of sample means (sd(means))
# = se.s

# when iter value goes to 
# infinite value:
# mean(means) = mean(p1) 
# and
# sd(means) = sd(p1) / sqrt(s.size)
# that is, se.s = se.z
# This is called CLT (Central Limit Theorem)
# see http://commres.net/wiki/cetral_limit_theorem

mean(means)
mean(p1)
sd(means)
var(p1) 
# remember we started talking sample size 10
sqrt(var(p1)/s.size)
se.z

sd(means)
se.s 
se.z

# because CLT
loz1 <- mean(p1)-se.z
hiz1 <- mean(p1)+se.z
loz2 <- mean(p1)-2*se.z
hiz2 <- mean(p1)+2*se.z
loz3 <- mean(p1)-3*se.z
hiz3 <- mean(p1)+3*se.z

c(lo1, loz1)
c(lo2, loz2)
c(lo3, loz3)

c(hi1, hiz1)
c(hi2, hiz2)
c(hi3, hiz3)


hist(means, breaks=50,
     xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)), 
     col = rgb(1, 1, 1, .5))
abline(v=mean(means), col="black", lwd=3)
# abline(v=mean(p2), colo='darkgreen', lwd=3)
abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3), 
       col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"), 
       lwd=2)

round(c(lo1, hi1))
round(c(lo2, hi2))
round(c(lo3, hi3))

round(c(loz1, hiz1))
round(c(loz2, hiz2))
round(c(loz3, hiz3))

m.sample.i.got <- mean(means)+ 1.5*sd(means)
m.sample.i.got

hist(means, breaks=30, 
     xlim = c(mean(means)-7*sd(means), mean(means)+10*sd(means)), 
     col = rgb(1, 1, 1, .5))
abline(v=mean(means), col="black", lwd=3)
abline(v=m.sample.i.got, col='darkgreen', lwd=3)

# what is the probablity of getting
# greater than 
# m.sample.i.got?
m.sample.i.got
pnorm(m.sample.i.got, mean(means), sd(means), lower.tail = F)
pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F)

# then, what is the probabilty of getting 
# greater than m.sample.i.got and
# less than corresponding value, which is
# mean(means) - m.sample.i.got - mean(means)
# (green line)
tmp <- mean(means) - (m.sample.i.got - mean(means))
abline(v=tmp, col='red', lwd=3)
2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
m.sample.i.got

### one more time 
# this time, with a story
mean(p2)
sd(p2)
s.from.p2 <- sample(p2, s.size)
m.s.from.p2 <- mean(s.from.p2)
m.s.from.p2

se.s
se.z
sd(means)

tmp <- mean(means) - (m.s.from.p2 - mean(means))
tmp 

hist(means, breaks=30, 
     xlim = c(tmp-4*sd(means), m.s.from.p2+4*sd(means)), 
     col = rgb(1, 1, 1, .5))
abline(v=mean(means), col="black", lwd=3)
abline(v=m.s.from.p2, col='blue', lwd=3)

# what is the probablity of getting
# greater than 
# m.sample.i.got?
m.s.from.p2
pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)

# then, what is the probabilty of getting 
# greater than m.sample.i.got and
# less than corresponding value, which is
# mean(means) - m.sample.i.got - mean(means)
# (green line)
abline(v=tmp, col='red', lwd=3)
2 * pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)

se.z
sd(s.from.p2)/sqrt(s.size)
se.z.adjusted <- sqrt(var(s.from.p2)/s.size)
se.z.adjusted
2 * pnorm(m.s.from.p2, mean(p1), se.z.adjusted, 
          lower.tail = F)

z.cal <- (m.s.from.p2 - mean(p1))/se.z.adjusted
z.cal
pnorm(z.cal, lower.tail = F)*2


pt(z.cal, 49, lower.tail = F)*2
t.test(s.from.p2, mu=mean(p1), var.equal = T)
</code>

====== output ======
<code>
> rm(list=ls())
> rnorm2 <- function(n,mean,sd){ 
+   mean+sd*scale(rnorm(n)) 
+ }
> 
> se <- function(sample) {
+   sd(sample)/sqrt(length(sample))
+ }
> 
> ss <- function(x) {
+   sum((x-mean(x))^2)
+ }
> 
> ################################
> N.p <- 1000000
> m.p <- 100
> sd.p <- 10
> 
> p1 <- rnorm2(N.p, m.p, sd.p)
> mean(p1)
[1] 100
> sd(p1)
[1] 10
> 
> p2 <- rnorm2(N.p, m.p+5, sd.p)
> mean(p2)
[1] 105
> sd(p2)
[1] 10
> 
> m.p1 <- mean(p1)
> sd.p1 <- sd(p1)
> var(p1)
     [,1]
[1,]  100
> 
> 
> hist(p1, breaks=50, col = rgb(1, 1, 1, 0.5),
+      main = "histogram of p1 and p2",)
> abline(v=mean(p1), col="black", lwd=3)
> hist(p2, add=T, breaks=50, col=rgb(1,1,.5,.5))
> abline(v=mean(p2), col="red", lwd=3)
> 
> 
> hist(p1, breaks=50, col=rgb(0,.5,.5,.5))
> abline(v=mean(p1),lwd=2)
> abline(v=mean(p1)-sd(p1), lwd=2)
> abline(v=mean(p1)+sd(p1), lwd=2)
> abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
> abline(v=c(m.p1-3*sd.p1, m.p1+3*sd.p1), lwd=2, col='green')
> 
> # area bet black = 68%
> # between red = 95%
> # between green = 99%
> 
> pnorm(m.p1+sd.p1, m.p1, sd.p1)
[1] 0.8413447
> pnorm(m.p1-sd.p1, m.p1, sd.p1)
[1] 0.1586553
> pnorm(m.p1+sd.p1, m.p1, sd.p1) - 
+   pnorm(m.p1-sd.p1, m.p1, sd.p1)
[1] 0.6826895
> 
> pnorm(m.p1+2*sd.p1, m.p1, sd.p1) - 
+   pnorm(m.p1-2*sd.p1, m.p1, sd.p1)
[1] 0.9544997
> 
> pnorm(m.p1+3*sd.p1, m.p1, sd.p1) - 
+   pnorm(m.p1-3*sd.p1, m.p1, sd.p1)
[1] 0.9973002
> 
> m.p1
[1] 100
> sd.p1
[1] 10
> 
> (m.p1-m.p1)/sd.p1
[1] 0
> ((m.p1-sd.p1) - m.p1) / sd.p1
[1] -1
> (120-100)/10 
[1] 2
> pnorm(1)-pnorm(-1)
[1] 0.6826895
> pnorm(2)-pnorm(-2)
[1] 0.9544997
> pnorm(3)-pnorm(3)
[1] 0
> 
> 1-pnorm(-2)*2
[1] 0.9544997
> pnorm(2)-pnorm(-2)
[1] 0.9544997
> 
> pnorm(120, 100, 10)
[1] 0.9772499
> pnorm(2)-pnorm(-2)
[1] 0.9544997
> 
> zscore <- (120-100)/10
> pnorm(zscore)-pnorm(-zscore)
[1] 0.9544997
> zscore
[1] 2
> 
> # reminder. 
> pnorm(-1)
[1] 0.1586553
> pnorm(1, 0, 1, lower.tail = F)
[1] 0.1586553
> pnorm(110, 100, 10, lower.tail = F)
[1] 0.1586553
> zscore <- (110-100)/10
> pnorm(zscore, lower.tail = F)
[1] 0.1586553
> 
> pnorm(118, 100, 10, lower.tail = F)
[1] 0.03593032
> pnorm(18/10, lower.tail = F)
[1] 0.03593032
> 
> z.p1 <- (p1-mean(p1))/sd(p1)
> mean(z.p1)
[1] 8.189457e-18
> round(mean(z.p1),10)
[1] 0
> sd(z.p1)
[1] 1
> pnorm(1.8)-pnorm(-1.8)
[1] 0.9281394
> 
> hist(z.p1, breaks=50, col=rgb(1,0,0,0))
> abline(v=c(m.p1, -1.8, 1.8), col='red')
> 1-(pnorm(1.8)-pnorm(-1.8))
[1] 0.07186064
> 
> 
> pnorm(1)-pnorm(-1)
[1] 0.6826895
> 1-(pnorm(-1)*2)
[1] 0.6826895
> 
> pnorm(2)-pnorm(-2)
[1] 0.9544997
> 1-(pnorm(-2)*2)
[1] 0.9544997
> 
> 1-(pnorm(-3)*2)
[1] 0.9973002
> 
> #
> hist(p1, breaks=50, col=rgb(.9,.9,.9,.9))
> abline(v=mean(p1),lwd=2)
> abline(v=mean(p1)-sd(p1), lwd=2)
> abline(v=mean(p1)+sd(p1), lwd=2)
> abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
> abline(v=c(m.p1-3*sd.p1, m.p1+3*sd.p1), lwd=2, col='green')
> 
> # 68%
> a <- qnorm(.32/2)
> b <- qnorm(1-.32/2)
> c(a, b)
[1] -0.9944579  0.9944579
> c(-1, 1)
[1] -1  1
> # note that
> .32/2
[1] 0.16
> pnorm(-1)
[1] 0.1586553
> qnorm(.32/2)
[1] -0.9944579
> qnorm(pnorm(-1))
[1] -1
> 
> # 95%
> c <- qnorm(.05/2)
> d <- qnorm(1-.05/2)
> c(c, d)
[1] -1.959964  1.959964
> c(-2,2)
[1] -2  2
> 
> # 99% 
> e <- qnorm(.01/2)
> f <- qnorm(1-.01/2)
> c(e,f)
[1] -2.575829  2.575829
> c(-3,3)
[1] -3  3
> 
> 
> pnorm(b)-pnorm(a)
[1] 0.68
> c(a, b)
[1] -0.9944579  0.9944579
> pnorm(d)-pnorm(c)
[1] 0.95
> c(c, d)
[1] -1.959964  1.959964
> pnorm(f)-pnorm(e)
[1] 0.99
> c(e, f)
[1] -2.575829  2.575829
> 
> qnorm(.5)
[1] 0
> qnorm(1)
[1] Inf
> 
> 
> ################################
> s.size <- 50
> 
> means.temp <- c()
> s1 <- sample(p1, s.size, replace = T)
> mean(s1)
[1] 98.76098
> means.temp <- append(means.temp, mean(s1))
> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
> means.temp
[1]  98.76098  99.90935  99.29643  99.66014 101.93822
> 
> iter <- 1000000
> # means <- c()
> means <- rep(NA, iter)
> for (i in 1:iter) {
+   # means <- append(means, mean(sample(p1, s.size, replace = T)))
+   means[i] <- mean(sample(p1, s.size, replace = T))
+ }
> length(means)
[1] 1000000
> mean(means)
[1] 99.99824
> sd(means)
[1] 1.414035
> se.s <- sd(means)
> 
> hist(means, breaks=50, 
+      xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)), 
+      col=rgb(1, 1, 1, .5))
> abline(v=mean(means), col="black", lwd=3)
> # now we want to get sd of this distribution
> lo1 <- mean(means)-se.s
> hi1 <- mean(means)+se.s
> lo2 <- mean(means)-2*se.s
> hi2 <- mean(means)+2*se.s
> lo3 <- mean(means)-3*se.s
> hi3 <- mean(means)+3*se.s
> abline(v=mean(means), col="black", lwd=2)
> # abline(v=mean(p2), colo='darkgreen', lwd=2)
> abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3), 
+        col=c("red","red", "blue", "blue", "orange", "orange"), 
+        lwd=2)
> 
> # meanwhile . . . .
> se.s
[1] 1.414035
> 
> se.z <- sqrt(var(p1)/s.size)
> se.z <- c(se.z)
> se.z
[1] 1.414214
> 
> se.z.adjusted <- sqrt(var(s1)/s.size)
> se.z.adjusted
[1] 1.370464
> 
> # sd of sample means (sd(means))
> # = se.s
> 
> # when iter value goes to 
> # infinite value:
> # mean(means) = mean(p1) 
> # and
> # sd(means) = sd(p1) / sqrt(s.size)
> # that is, se.s = se.z
> # This is called CLT (Central Limit Theorem)
> # see http://commres.net/wiki/cetral_limit_theorem
> 
> mean(means)
[1] 99.99824
> mean(p1)
[1] 100
> sd(means)
[1] 1.414035
> var(p1) 
     [,1]
[1,]  100
> # remember we started talking sample size 10
> sqrt(var(p1)/s.size)
         [,1]
[1,] 1.414214
> se.z
[1] 1.414214
> 
> sd(means)
[1] 1.414035
> se.s 
[1] 1.414035
> se.z
[1] 1.414214
> 
> # because CLT
> loz1 <- mean(p1)-se.z
> hiz1 <- mean(p1)+se.z
> loz2 <- mean(p1)-2*se.z
> hiz2 <- mean(p1)+2*se.z
> loz3 <- mean(p1)-3*se.z
> hiz3 <- mean(p1)+3*se.z
> 
> c(lo1, loz1)
[1] 98.58421 98.58579
> c(lo2, loz2)
[1] 97.17017 97.17157
> c(lo3, loz3)
[1] 95.75614 95.75736
> 
> c(hi1, hiz1)
[1] 101.4123 101.4142
> c(hi2, hiz2)
[1] 102.8263 102.8284
> c(hi3, hiz3)
[1] 104.2403 104.2426
> 
> 
> hist(means, breaks=50,
+      xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)), 
+      col = rgb(1, 1, 1, .5))
> abline(v=mean(means), col="black", lwd=3)
> # abline(v=mean(p2), colo='darkgreen', lwd=3)
> abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3), 
+        col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"), 
+        lwd=2)
> 
> round(c(lo1, hi1))
[1]  99 101
> round(c(lo2, hi2))
[1]  97 103
> round(c(lo3, hi3))
[1]  96 104
> 
> round(c(loz1, hiz1))
[1]  99 101
> round(c(loz2, hiz2))
[1]  97 103
> round(c(loz3, hiz3))
[1]  96 104
> 
> m.sample.i.got <- mean(means)+ 1.5*sd(means)
> m.sample.i.got
[1] 102.1193
> 
> hist(means, breaks=30, 
+      xlim = c(mean(means)-7*sd(means), mean(means)+10*sd(means)), 
+      col = rgb(1, 1, 1, .5))
> abline(v=mean(means), col="black", lwd=3)
> abline(v=m.sample.i.got, col='darkgreen', lwd=3)
> 
> # what is the probablity of getting
> # greater than 
> # m.sample.i.got?
> m.sample.i.got
[1] 102.1193
> pnorm(m.sample.i.got, mean(means), sd(means), lower.tail = F)
[1] 0.0668072
> pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F)
[1] 0.0669929
> 
> # then, what is the probabilty of getting 
> # greater than m.sample.i.got and
> # less than corresponding value, which is
> # mean(means) - m.sample.i.got - mean(means)
> # (green line)
> tmp <- mean(means) - (m.sample.i.got - mean(means))
> abline(v=tmp, col='red', lwd=3)
> 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
[1] 0.1339368
> m.sample.i.got
[1] 102.1193
> 
> ### one more time 
> # this time, with a story
> mean(p2)
[1] 105
> sd(p2)
[1] 10
> s.from.p2 <- sample(p2, s.size)
> m.s.from.p2 <- mean(s.from.p2)
> m.s.from.p2
[1] 103.1283
> 
> se.s
[1] 1.414035
> se.z
[1] 1.414214
> sd(means)
[1] 1.414035
> 
> tmp <- mean(means) - (m.s.from.p2 - mean(means))
> tmp 
[1] 96.86822
> 
> hist(means, breaks=30, 
+      xlim = c(tmp-4*sd(means), m.s.from.p2+4*sd(means)), 
+      col = rgb(1, 1, 1, .5))
> abline(v=mean(means), col="black", lwd=3)
> abline(v=m.s.from.p2, col='blue', lwd=3)
> 
> # what is the probablity of getting
> # greater than 
> # m.sample.i.got?
> m.s.from.p2
[1] 103.1283
> pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)
[1] 0.01348266
> 
> # then, what is the probabilty of getting 
> # greater than m.sample.i.got and
> # less than corresponding value, which is
> # mean(means) - m.sample.i.got - mean(means)
> # (green line)
> abline(v=tmp, col='red', lwd=3)
> 2 * pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)
[1] 0.02696533
> 
> se.z
[1] 1.414214
> sd(s.from.p2)/sqrt(s.size)
[1] 1.414296
> se.z.adjusted <- sqrt(var(s.from.p2)/s.size)
> se.z.adjusted
[1] 1.414296
> 2 * pnorm(m.s.from.p2, mean(p1), se.z.adjusted, 
+           lower.tail = F)
[1] 0.02697421
> 
> z.cal <- (m.s.from.p2 - mean(p1))/se.z.adjusted
> z.cal
[1] 2.211891
> pnorm(z.cal, lower.tail = F)*2
[1] 0.02697421
> 
> 
> pt(z.cal, 49, lower.tail = F)*2
[1] 0.03166797
> t.test(s.from.p2, mu=mean(p1), var.equal = T)

	One Sample t-test

data:  s.from.p2
t = 2.2119, df = 49, p-value = 0.03167
alternative hypothesis: true mean is not equal to 100
95 percent confidence interval:
 100.2861 105.9704
sample estimates:
mean of x 
 103.1283 

>
</code>

====== T-test sum up ======
<code>

rm(list=ls())
rnorm2 <- function(n,mean,sd){ 
  mean+sd*scale(rnorm(n)) 
}
se <- function(sample) {
  sd(sample)/sqrt(length(sample))
}
ss <- function(x) {
  sum((x-mean(x))^2)
}

N.p <- 1000000
m.p <- 100
sd.p <- 10


p1 <- rnorm2(N.p, m.p, sd.p)
mean(p1)
sd(p1)

p2 <- rnorm2(N.p, m.p+10, sd.p)
mean(p2)
sd(p2)

hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
     main = "histogram of p1 and p2",)
abline(v=mean(p1), col="black", lwd=3)
hist(p2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
abline(v=mean(p2), col="red", lwd=3)

s.size <- 1000
s2 <- sample(p2, s.size)
mean(s2)
sd(s2)

se.z <- sqrt(var(p1)/s.size)
se.z <- c(se.z)
se.z.range <- c(-2*se.z,2*se.z)
se.z.range

mean(p1)+se.z.range
mean(s2)

z.cal <- (mean(s2) - mean(p1)) / se.z
z.cal
pnorm(z.cal, lower.tail = F) * 2

z.cal 

# principles . . . 
# distribution of sample means 
iter <- 100000
means <- c()
for (i in 1:iter) {
  m.of.s <- mean(sample(p1, s.size))
  means <- append(means, m.of.s)
}

hist(means, 
     xlim = c(mean(means)-3*sd(means), mean(s2)+5), 
     col = rgb(1, 0, 0, 0.5))
abline(v=mean(p1), col="black", lwd=3)
abline(v=mean(s2), 
       col="blue", lwd=3)
lo1 <- mean(p1)-se.z
hi1 <- mean(p1)+se.z
lo2 <- mean(p1)-2*se.z
hi2 <- mean(p1)+2*se.z
abline(v=c(lo1, hi1, lo2, hi2), 
       col=c("green","green", "brown", "brown"), 
       lwd=2)
se.z
c(lo2, hi2)
pnorm(z.cal, lower.tail = F) * 2


# Note that we use sqrt(var(s2)/s.size)
# as our se value instread of 
# sqrt(var(p1)/s.size)
# This is a common practice for R
# In fact, some z.test (made by someone)
# function uses the former rather than
# latter.

sqrt(var(p1)/s.size)
se.z

sqrt(var(s2)/s.size)
se(s2)

t.cal.os <- (mean(s2) - mean(p1))/ se(s2)
z.cal <- (mean(s2) - mean(p1))/ se.z
t.cal.os
z.cal 

df.s2 <- length(s2)-1
df.s2
p.t.os <- pt(abs(t.cal.os), df.s2, lower.tail = F) * 2
p.t.os
t.out <- t.test(s2, mu=mean(p1))

library(BSDA)
z.out <- z.test(s2, p1, sigma.x = sd(s2), sigma.y = sd(p1))

z.out$statistic # se.z 대신에 se(s2) 값으로 구한 z 값
z.cal # se.z으로 (sqrt(var(p1)/s.size)값) 구한 z 값

t.out$statistic # se(s2)를 분모로 하여 구한 t 값
t.cal.os # se(s2)를 이용하여 손으로 구한 t 값

# But, after all, we use t.test method regardless of 
# variation

hist(means, 
     xlim = c(mean(means)-3*sd(means), mean(s2)+5), 
     col = rgb(1, 0, 0, 0.5))
abline(v=mean(p1), col="black", lwd=3)
abline(v=mean(s2), 
       col="blue", lwd=3)
lo1 <- mean(p1)-se.z
hi1 <- mean(p1)+se.z
lo2 <- mean(p1)-2*se.z
hi2 <- mean(p1)+2*se.z
abline(v=c(lo1, hi1, lo2, hi2), 
       col=c("green","green", "brown", "brown"), 
       lwd=2)

# difference between black and blue line
# divided by
# se(s2) (= random difference)
# t.value
mean(s2)
mean(p1)
diff <- mean(s2)-mean(p1)
diff
se(s2)
diff/se(s2)

t.cal.os

########################
# 2 sample t-test
########################
# 가정. 아래에서 추출하는 두 
# 샘플의 모집단의 파라미터를
# 모른다. 
s1 <- sample(p1, s.size)
s2 <- sample(p2, s.size)

mean(s1)
mean(s2)
ss(s1)
ss(s2)
df.s1 <- length(s1)-1
df.s2 <- length(s2)-1
df.s1
df.s2

pooled.variance <- (ss(s1)+ss(s2))/(df.s1+df.s2)
pooled.variance
se.ts <- sqrt((pooled.variance/length(s1))+(pooled.variance/length(s2)))
se.ts
t.cal.ts <- (mean(s1)-mean(s2))/se.ts
t.cal.ts
p.val.ts <- pt(abs(t.cal.ts), df=df.s1+df.s2, lower.tail = F) * 2
p.val.ts

t.test(s1, s2, var.equal = T)

se(s1)
se(s2)

mean(s1)+c(-se(s1)*2, se(s1)*2)
mean(s2)+c(-se(s2)*2, se(s2)*2)

mean(p1)
mean(p2)

hist(s1, breaks=50,
     col = rgb(1, 0, 0, 0.5))
hist(s2, breaks=50, add=T, col=rgb(0,0,1,1))
abline(v=mean(s1), col="green", lwd=3)
# hist(s2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
abline(v=mean(s2), col="lightblue", lwd=3)

diff <- mean(s1)-mean(s2)
se.ts
diff/se.ts

####
# repeated measure t-test
# we can use the above case 
# pop paramter unknown
# two consecutive measurement 
# for the same sample 

t1 <- s1
t2 <- s2
mean(t1)
mean(t2)
diff.s <- t1 - t2
diff.s
t.cal.rm <- mean(diff.s)/se(diff.s)
t.cal.rm
p.val.rm <- pt(abs(t.cal.rm), length(s1)-1, lower.tail = F) * 2
p.val.rm
t.test(s1, s2, paired = T)

# create multiple histogram
s.all <- c(s1,s2)
mean(s.all)
hist(s1, col='grey', breaks=50, xlim=c(50, 150))
hist(s2, col='darkgreen', breaks=50, add=TRUE)
abline(v=c(mean(s.all)), 
       col=c("red"), lwd=3)
abline(v=c(mean(s1), mean(s2)), 
       col=c("black", "green"), lwd=3)

comb <- data.frame(s1,s2)
dat <- stack(comb)
head(dat)

m.tot <- mean(s.all)
m.s1 <- mean(s1)
m.s2 <- mean(s2)

ss.tot <- ss(s.all)
ss.s1 <- ss(s1)
ss.s2 <- ss(s2)

df.tot <- length(s.all)-1
df.s1 <- length(s1)-1
df.s2 <- length(s2)-1

ms.tot <- var(s.all)
ms.tot
ss.tot/df.tot

var(s1)
ss.s1 / df.s1

var(s2)
ss.s2 / df.s2

ss.b.s1 <- length(s1) * ((m.tot - m.s1)^2)
ss.b.s2 <- length(s2) * ((m.tot - m.s1)^2)
ss.bet <- ss.b.s1+ss.b.s2
ss.bet

ss.wit <- ss.s1 + ss.s2
ss.wit

ss.bet + ss.wit
ss.tot

library(dplyr)
# df.bet <- length(unique(dat)) - 1
df.bet <- nlevels(dat$ind) - 1
df.wit <- df.s1+df.s2
df.bet
df.wit
df.bet+df.wit
df.tot 

ms.bet <- ss.bet / df.bet
ms.wit <- ss.wit / df.wit
ms.bet
ms.wit

f.cal <- ms.bet / ms.wit
f.cal
pf(f.cal, df.bet, df.wit, lower.tail = F)


f.test <-  aov(dat$values~ dat$ind, data = dat)
summary(f.test)

sqrt(f.cal)
t.cal.ts

# this is anova after all. 
</code>

====== output ======
<code>
> rm(list=ls())
> rnorm2 <- function(n,mean,sd){ 
+   mean+sd*scale(rnorm(n)) 
+ }
> se <- function(sample) {
+   sd(sample)/sqrt(length(sample))
+ }
> ss <- function(x) {
+   sum((x-mean(x))^2)
+ }
> 
> N.p <- 1000000
> m.p <- 100
> sd.p <- 10
> 
> 
> p1 <- rnorm2(N.p, m.p, sd.p)
> mean(p1)
[1] 100
> sd(p1)
[1] 10
> 
> p2 <- rnorm2(N.p, m.p+10, sd.p)
> mean(p2)
[1] 110
> sd(p2)
[1] 10
> 
> hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
+      main = "histogram of p1 and p2",)
> abline(v=mean(p1), col="black", lwd=3)
> hist(p2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
> abline(v=mean(p2), col="red", lwd=3)
> 
> s.size <- 1000
> s2 <- sample(p2, s.size)
> mean(s2)
[1] 110.4892
> sd(s2)
[1] 10.36614
> 
> se.z <- sqrt(var(p1)/s.size)
> se.z <- c(se.z)
> se.z.range <- c(-2*se.z,2*se.z)
> se.z.range
[1] -0.6324555  0.6324555
> 
> mean(p1)+se.z.range
[1]  99.36754 100.63246
> mean(s2)
[1] 110.4892
> 
> z.cal <- (mean(s2) - mean(p1)) / se.z
> z.cal
[1] 33.16976
> pnorm(z.cal, lower.tail = F) * 2
[1] 2.93954e-241
> 
> z.cal 
[1] 33.16976
> 
> # principles . . . 
> # distribution of sample means 
> iter <- 100000
> means <- c()
> for (i in 1:iter) {
+   m.of.s <- mean(sample(p1, s.size))
+   means <- append(means, m.of.s)
+ }
> 
> hist(means, 
+      xlim = c(mean(means)-3*sd(means), mean(s2)+5), 
+      col = rgb(1, 0, 0, 0.5))
> abline(v=mean(p1), col="black", lwd=3)
> abline(v=mean(s2), 
+        col="blue", lwd=3)
> lo1 <- mean(p1)-se.z
> hi1 <- mean(p1)+se.z
> lo2 <- mean(p1)-2*se.z
> hi2 <- mean(p1)+2*se.z
> abline(v=c(lo1, hi1, lo2, hi2), 
+        col=c("green","green", "brown", "brown"), 
+        lwd=2)
> se.z
[1] 0.3162278
> c(lo2, hi2)
[1]  99.36754 100.63246
> pnorm(z.cal, lower.tail = F) * 2
[1] 2.93954e-241
> 
> 
> # Note that we use sqrt(var(s2)/s.size)
> # as our se value instread of 
> # sqrt(var(p1)/s.size)
> # This is a common practice for R
> # In fact, some z.test (made by someone)
> # function uses the former rather than
> # latter.
> 
> sqrt(var(p1)/s.size)
          [,1]
[1,] 0.3162278
> se.z
[1] 0.3162278
> 
> sqrt(var(s2)/s.size)
[1] 0.3278061
> se(s2)
[1] 0.3278061
> 
> t.cal.os <- (mean(s2) - mean(p1))/ se(s2)
> z.cal <- (mean(s2) - mean(p1))/ se.z
> t.cal.os
[1] 31.99818
> z.cal 
[1] 33.16976
> 
> df.s2 <- length(s2)-1
> df.s2
[1] 999
> p.t.os <- pt(abs(t.cal.os), df.s2, lower.tail = F) * 2
> p.t.os
[1] 3.162139e-155
> t.out <- t.test(s2, mu=mean(p1))
> 
> library(BSDA)
필요한 패키지를 로딩중입니다: lattice

다음의 패키지를 부착합니다: ‘BSDA’

The following object is masked from ‘package:datasets’:

    Orange

경고메시지(들): 
패키지 ‘BSDA’는 R 버전 4.3.3에서 작성되었습니다 
> z.out <- z.test(s2, p1, sigma.x = sd(s2), sigma.y = sd(p1))
> 
> z.out$statistic # se.z 대신에 se(s2) 값으로 구한 z 값
      z 
31.9833 
> z.cal # se.z으로 (sqrt(var(p1)/s.size)값) 구한 z 값
[1] 33.16976
> 
> t.out$statistic # se(s2)를 분모로 하여 구한 t 값
       t 
31.99818 
> t.cal.os # se(s2)를 이용하여 손으로 구한 t 값
[1] 31.99818
> 
> # But, after all, we use t.test method regardless of 
> # variation
> 
> hist(means, 
+      xlim = c(mean(means)-3*sd(means), mean(s2)+5), 
+      col = rgb(1, 0, 0, 0.5))
> abline(v=mean(p1), col="black", lwd=3)
> abline(v=mean(s2), 
+        col="blue", lwd=3)
> lo1 <- mean(p1)-se.z
> hi1 <- mean(p1)+se.z
> lo2 <- mean(p1)-2*se.z
> hi2 <- mean(p1)+2*se.z
> abline(v=c(lo1, hi1, lo2, hi2), 
+        col=c("green","green", "brown", "brown"), 
+        lwd=2)
> 
> # difference between black and blue line
> # divided by
> # se(s2) (= random difference)
> # t.value
> mean(s2)
[1] 110.4892
> mean(p1)
[1] 100
> diff <- mean(s2)-mean(p1)
> diff
[1] 10.4892
> se(s2)
[1] 0.3278061
> diff/se(s2)
[1] 31.99818
> 
> t.cal.os
[1] 31.99818
> 
> ########################
> # 2 sample t-test
> ########################
> # 가정. 아래에서 추출하는 두 
> # 샘플의 모집단의 파라미터를
> # 모른다. 
> s1 <- sample(p1, s.size)
> s2 <- sample(p2, s.size)
> 
> mean(s1)
[1] 100.7138
> mean(s2)
[1] 109.8804
> ss(s1)
[1] 108288.9
> ss(s2)
[1] 100751
> df.s1 <- length(s1)-1
> df.s2 <- length(s2)-1
> df.s1
[1] 999
> df.s2
[1] 999
> 
> pooled.variance <- (ss(s1)+ss(s2))/(df.s1+df.s2)
> pooled.variance
[1] 104.6246
> se.ts <- sqrt((pooled.variance/length(s1))+(pooled.variance/length(s2)))
> se.ts
[1] 0.4574376
> t.cal.ts <- (mean(s1)-mean(s2))/se.ts
> t.cal.ts
[1] -20.03892
> p.val.ts <- pt(abs(t.cal.ts), df=df.s1+df.s2, lower.tail = F) * 2
> p.val.ts
[1] 1.522061e-81
> 
> t.test(s1, s2, var.equal = T)

	Two Sample t-test

data:  s1 and s2
t = -20.039, df = 1998, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10.06366  -8.26945
sample estimates:
mean of x mean of y 
 100.7138  109.8804 

> 
> se(s1)
[1] 0.3292374
> se(s2)
[1] 0.3175719
> 
> mean(s1)+c(-se(s1)*2, se(s1)*2)
[1] 100.0553 101.3723
> mean(s2)+c(-se(s2)*2, se(s2)*2)
[1] 109.2452 110.5155
> 
> mean(p1)
[1] 100
> mean(p2)
[1] 110
> 
> hist(s1, breaks=50,
+      col = rgb(1, 0, 0, 0.5))
> hist(s2, breaks=50, add=T, col=rgb(0,0,1,1))
> abline(v=mean(s1), col="green", lwd=3)
> # hist(s2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
> abline(v=mean(s2), col="lightblue", lwd=3)
> 
> diff <- mean(s1)-mean(s2)
> se.ts
[1] 0.4574376
> diff/se.ts
[1] -20.03892
> 
> ####
> # repeated measure t-test
> # we can use the above case 
> # pop paramter unknown
> # two consecutive measurement 
> # for the same sample 
> 
> t1 <- s1
> t2 <- s2
> mean(t1)
[1] 100.7138
> mean(t2)
[1] 109.8804
> diff.s <- t1 - t2
> diff.s
   [1] -11.22445947  -3.01532631  -3.47460077 -14.77286374  -9.74734969   2.33544374  -1.53761723
   [8]   5.16258644 -40.75608052  -4.94498515  13.05099703  -0.24847628  -6.26581498  -8.14116058
  [15]   5.49152742   6.15799638  -4.65589524  -7.93104859 -15.28488198  -5.99182327 -15.35591992
  [22]  -1.00704696  18.27285355 -12.09936359  -4.72393289  -4.98301049 -20.95452153   1.62363928
  [29] -10.58244069 -21.50608157 -53.61230898  -3.85198479 -42.61929736  -6.80266370 -22.92704580
  [36]   3.01745740 -19.37131113 -27.82551902 -10.05485425 -25.12701225 -12.93162558  -7.55706006
  [43] -19.16855657  -4.81878797   1.64397602 -28.64658004  16.36241227   8.73170802 -11.56090742
  [50]   3.21799642 -39.37233043  -9.96051946 -19.11232333 -34.53077051  -4.85780005  -9.52501389
  [57]  -8.28743679 -38.33995044 -50.60884456  -3.43450084  -0.85381393 -13.30971467  10.13049966
  [64]   8.65616917 -29.75453733 -25.40674843 -24.98197786 -12.92901371  15.80168803  -8.67599446
  [71] -20.50728324 -19.37275012 -23.27866089 -11.74962053 -34.35317908 -26.10460199   8.59009957
  [78] -24.79252799   3.09475727 -19.13505970 -11.72561867 -33.79775614  -6.00167910 -25.03263480
  [85] -23.66447959  -8.54416282  -6.89905337 -10.45234583 -34.67182776 -37.20205500   6.60270378
  [92] -21.22842221 -11.68774036  -8.71535203 -15.55746542   2.88009050  -5.51543509   1.21606420
  [99]  19.63435733 -14.40578016 -11.24172079  10.21723258 -23.41564885  27.60247565  -8.28684078
 [106]  17.72472594  -0.29977586 -14.84142327 -20.25713391 -37.99518419 -27.68545647 -19.78976153
 [113] -10.23092427   0.40875267 -17.36077213 -24.53979674 -24.18070810 -24.13321556 -17.36615616
 [120] -31.50478963  -2.47101725 -22.14003910 -33.63875270 -19.91485505  -3.64251563 -30.62407901
 [127]  -7.91406849  -6.25389594   4.40651820 -19.08031290 -14.01489366 -31.30542657 -27.05443597
 [134]  -6.60642443   1.29892762 -11.73908399   1.96878666  -7.53666283 -42.78007247  -9.04715952
 [141]  14.05326537  -6.85724091 -16.02305236 -24.38581030  -9.91759245   4.75488243  -2.83181250
 [148]  -5.38371023 -19.62202451  -1.55000290 -14.86541899   2.95279749   3.82076747  -3.38868029
 [155]  -0.65101074  18.42861149 -20.55548459   3.02438240  21.23539589   3.32810800  -9.66847192
 [162]   2.48687983  -5.40571073  -1.53265391 -12.93542011  19.87564176 -10.69228781  12.80629134
 [169]   6.51132022 -16.96586244  10.88690774  -0.37382590 -14.69255590 -26.66941722  -3.67854181
 [176]   6.29443209 -10.77182585 -25.65187453   0.09825251   9.03908176  10.25414786  -0.45340800
 [183] -38.63379525 -16.58800530 -17.62672134   5.43887886  -4.95532070  -6.46372777  -3.14562140
 [190] -22.25276559   2.05941880 -44.33676979  -3.34190343 -19.04858920  -8.21990394 -25.53625536
 [197] -17.21120659  28.96404607   6.25767994  -6.17593254 -34.33503461   6.65350479  11.42897662
 [204]  -0.83715976 -28.46397824 -40.67262397  -6.01225907 -11.22598108 -10.92756008 -15.45671946
 [211]  -4.57060131 -27.19860432   5.68618678 -27.70257611 -27.77374648  -5.93312400  -9.37871992
 [218] -24.41623403  16.94244832 -30.46760860  -4.91996788   2.89031604   5.27074167 -23.91666746
 [225] -13.27091592  -7.99640540  15.26148582 -26.01138488 -28.57927092 -17.29274303 -17.07704891
 [232] -11.52528966  -5.50387909  -0.66159232 -11.50347650 -19.90680762 -11.09595230 -11.02710712
 [239] -13.34969773 -37.98584006  -0.95289265 -15.00431567  -1.22592809   7.40922588  11.45790664
 [246] -24.07983488   4.54606079  -0.54863357   6.56528626  -4.04491250 -17.13525622 -22.85976576
 [253] -12.30101864   7.01445235  -6.77058075 -10.69023765 -14.21289974   0.68743488 -22.13964282
 [260]  -4.93155960  -5.32992121  -9.24699990 -34.21542676 -28.10074867 -14.64483350 -21.94636738
 [267]   0.57190289 -32.90838279 -23.39341251  -8.52122572   7.61461839 -12.60688433  -8.17161329
 [274]  -7.73981345  -4.86979671  14.97509924 -17.25493992 -48.85010339  10.16448581  -4.34608694
 [281]  -4.73924884  19.07076764  -5.23571728 -28.73076387  -1.01530521  -1.14387890  -6.96197277
 [288] -10.10160879 -11.59352210   5.83294359 -10.03967522  -9.59761019 -15.88867160  -7.13643475
 [295]  -7.29391701 -31.93027109 -15.56526408  13.22678162 -11.18996097 -25.00719650 -12.64524490
 [302] -35.18133234  13.41791211 -10.87228845 -31.95732546 -18.32496165 -17.76804212 -14.54640557
 [309]  21.57859526   8.22556833 -11.51111550  -7.19669828 -28.98417620  -8.16801696  -3.70794736
 [316]  -6.13751897   8.57292816   0.83289680 -11.25989874  11.25387495 -21.86409747   8.11413189
 [323] -31.81200194 -12.67414744  -7.24837917  -9.76383734  -9.95068555 -17.43835347 -21.88852882
 [330]  -7.57724957 -39.97109963 -18.22405067  -3.35304894  -2.01422861  12.65855868  19.68084692
 [337] -17.60471709   0.17064467  -2.41707219  -8.99719317  -6.91454803   6.10676201  -6.24933210
 [344] -10.52672545  -8.28580388 -21.74741007 -20.10217608 -27.73408273 -15.76758951 -16.05537057
 [351]   7.28801425  14.00094463 -32.21492728 -26.60734855   5.34256687  -2.73678515 -11.37626522
 [358] -23.18199896 -18.58488442 -26.09162223  -8.11379506 -16.24314767  -7.63202196  -7.70819353
 [365]  15.27344158  -6.89298171 -14.69008686 -10.61871285 -15.97716535   2.04244654 -20.72377059
 [372]  -3.33710996  -8.70719328   1.11210610  17.29240572   3.03254095  -0.17471362  -0.03964601
 [379] -37.25886118   1.80635404 -18.13062850 -12.08353080   3.17736630  11.60663212  14.40230421
 [386] -10.17057116 -11.45984282 -11.29130871 -15.60613249   3.94447987  21.39005539  38.66131508
 [393] -21.77681110  12.35832123  20.97222557 -30.89221977  26.68094540  13.28471640   5.74956285
 [400]  -7.60656480  17.47154758   2.14422376 -25.38609321  -7.75483605  12.58051687  -0.12179751
 [407] -17.52668230   1.92572455  -0.40532799 -11.06357792 -11.11874742 -10.35835894 -18.01080311
 [414] -14.30821541 -10.88427284   3.35560021 -13.46512417   6.04047559 -14.68666762   1.68951313
 [421]   4.31329242   7.50742041 -26.79115307 -32.01096580  28.21960400   3.60371811  14.05530287
 [428] -23.35919604   4.65486577 -42.05164227   5.90447867 -22.06140361 -36.42899364 -37.85288612
 [435]  -7.63484440  -7.99917501   5.73410720 -17.24124704  25.18142781 -22.27250215   3.31690400
 [442] -19.75974381  13.54264541 -18.86771849 -37.61540899   1.88862667 -10.85162959 -15.90454635
 [449]  -2.21038773  -1.64057323  -5.05984647 -11.33623396  -3.23993046  -7.07393112 -12.97757803
 [456]   2.81189108 -14.48709633 -26.49136729 -27.92148265  -3.32407602 -32.81165304  12.60241883
 [463] -19.77773032  -7.21862902 -30.05644966 -34.95205067  11.36207603  -1.88732702  -0.76534355
 [470]  -2.96997060 -26.06871359 -35.94190403  15.37335724 -29.72925340 -16.82830601 -32.14796579
 [477]  -7.49787668 -18.90903047  -4.95224428 -15.73171387  13.33385011   6.12192173 -12.11076372
 [484] -18.62856464 -25.58573774 -11.07888550 -16.77332405  -0.02979250   0.28045921   6.30825053
 [491]   2.66879530 -12.26308670  -3.34380103   0.57984817  -6.07862084  -0.47454107   0.18277721
 [498]  -1.49326079 -10.71424083 -18.55815468 -12.12523481 -10.32886876   1.21618496 -26.03417587
 [505]   4.21652961  -7.19317598  -6.49276222 -24.15078462 -16.62200664  18.65907236 -14.55468004
 [512]  12.03275099  13.63985493   1.94955005   4.24574372  -1.87352631  -9.90476053   4.92904187
 [519] -27.58434557  -3.22307117 -17.00433045 -24.92647628  -9.86853681  -3.12201969  -3.23781294
 [526] -44.00039083 -11.40638340 -10.42772214  11.56731373  34.18935566   0.22493781 -21.75192741
 [533] -27.11548294 -26.69534077  11.86361692 -18.87587506 -19.76960989  -6.72305933 -16.37844743
 [540]  28.42329924  -9.05042270  -5.94334753   0.80277394 -10.71835822  -7.39049563  -2.78087519
 [547] -18.87000360  -2.56690042  -2.31189557 -15.13438755   2.47712830  -7.13675100  -0.16789790
 [554] -17.69850124 -11.92887098  -7.41497303 -12.51272098   8.45361690 -30.63879338 -22.19654933
 [561]   0.37289814 -27.30901665 -42.04696582 -12.04972054 -20.74642541  -7.29054494 -11.52182026
 [568]   8.38372199 -18.57905859   7.62453900 -15.52892307  -1.83154260 -16.24286276  -5.74980635
 [575]   5.80843070   5.89082648 -12.38211486 -30.94945977 -34.27141760 -16.68862227 -12.74228148
 [582] -12.87197389 -27.14630137 -14.53291112  -8.24059195  -5.93523740   3.67320111 -14.73114416
 [589]   7.12333661 -20.54493654 -16.23259278  -7.23953837  13.66853283 -13.35719133  -3.19469244
 [596]  -8.63453073  -4.38937208  -4.25683530  -3.01229288  -9.46716386 -14.49889042  -4.80320886
 [603]  -5.93786510  -5.58582436 -38.52292640  -7.31871320 -10.10261534   3.55750183  -0.27859803
 [610]  13.47634673 -29.68104705  -4.28539812   3.67080445  -5.84424964  13.18400307  -7.89086592
 [617]  -1.33920025  -0.95308433 -29.91743460 -18.66903033 -20.15406028  -0.89348972 -12.41231804
 [624]   8.08680667 -34.26161156 -33.98947051  -7.80666121 -10.20912967 -28.95896039  -5.89738802
 [631]   2.32928389 -38.93521867 -15.56868160  -6.70412659  -2.57471692   2.46451660 -12.81558476
 [638]  19.76309298   2.12194484  24.84734206 -13.55620029 -17.87794609 -19.54477100 -19.73182713
 [645] -42.02480771 -13.42077241   9.19843679 -15.49829521 -35.16885995 -17.77559877  -2.28384147
 [652] -17.68303368   6.17812431   0.66249967   5.00542555  -2.26766815 -11.73064973  -6.75727090
 [659] -19.48957914 -28.88885682 -10.59583907 -32.20537872 -11.95661953 -23.77887711 -23.51551361
 [666] -17.32443690 -10.86310515  -8.51040349 -27.57610667 -26.85875525  -3.70329222 -20.50863173
 [673] -12.41671968  -5.15745402 -19.62849573   4.85082387 -29.88557391 -20.24914582 -27.73847105
 [680]  -0.12605203  -9.20839167 -51.56244236   3.17764075  13.96786787 -12.74398310   6.86534410
 [687] -21.75000477  14.10169236 -24.69667641 -20.26619614 -11.67028168  -5.14496708 -21.84000650
 [694] -34.30010114  -4.30214907 -15.81158253  -2.54412477   0.17601622 -22.22290730  -6.51460318
 [701] -19.46561809   9.62212347  -5.62354822  16.70312068  -7.16879691  -5.77420998  -2.36157455
 [708] -27.91638644 -12.61381331   7.45329002  -1.78749631  10.24888993  -2.76665687  -6.47189694
 [715] -20.55376627   1.03372077  -6.75380336 -21.29024889 -16.12342903 -36.81337018   1.75482644
 [722] -14.38944775  -4.27006397 -10.21581755   1.97016866   1.50969462  32.31451580   3.32233756
 [729]  -7.85868267   8.83356066  -3.54004596 -17.21481071 -29.58350979 -22.72248706 -27.86169027
 [736] -20.25705972  -1.67627671 -23.02237081 -20.13752529   6.07661361  -0.84839297  19.31619624
 [743] -13.32818441  -1.51206927 -16.05469364 -18.19320869  21.19248327 -26.85398142   2.82896396
 [750]   1.90853566 -10.76451371  -0.16368097  -5.02703204 -23.15483742  -5.12822113  -4.84245502
 [757] -19.16011286 -13.91801221 -11.66649472  15.04676653 -20.54651422  32.67150526 -11.37626079
 [764] -14.75337241  -0.46891630 -12.09854921 -10.75658868 -18.03655441  10.95871312  27.68695247
 [771]  -1.22676076 -15.78897397 -13.68374038  16.98138996 -15.57048042 -17.53983895 -32.33929466
 [778] -34.01869977  -2.46227514  -6.56500408 -17.04103052 -17.24440339  -6.68381805  -8.43674456
 [785]  -3.88372407 -29.28134174 -40.86613320 -18.64259952  -2.78880196  11.13938536  14.40875929
 [792]  -6.52858972 -38.92485161 -23.57441915 -25.59877817  19.51852353 -20.06650023  -4.32313864
 [799]   4.62056706  -5.18094527  -0.76429848  -2.98392902 -18.50300318 -29.63778693 -23.63235389
 [806]   5.68017122  14.33220812 -13.31317005 -10.81641168   5.22445430 -35.46250519   1.56327962
 [813] -14.60867384  10.29147319 -13.28593538   6.63825469 -14.52606348  17.45903705 -38.05095694
 [820] -13.93270232 -20.21993468  -1.96308711   1.80444271  12.16855255   9.52956342 -14.88194384
 [827] -29.04193544 -24.04102844 -14.01878071  -2.03269506   2.67151865  -5.20017572 -41.14943705
 [834]   8.96661691  28.12018815  -2.37196235 -13.46669223 -15.34687871 -19.92157033  16.10283716
 [841]   5.71060454  -1.87210810  17.82634786   6.46299261 -23.56325888  -9.31538158  13.08900119
 [848]   8.81863004 -16.87823373   2.57469446 -19.81326240  -1.08297141 -15.99656489 -12.78570251
 [855]  -8.53943328 -37.51286174  -5.80934175  29.94051347 -12.29916397  -0.84174744  -5.89053659
 [862] -30.93593593  -6.24638974   6.71567898   0.33777483  -6.43007412  -6.10032287 -18.90969351
 [869]   6.22885535   2.29565188 -25.72416278  -4.48305502  -5.77922453 -13.55585021 -23.84825362
 [876] -14.65449874 -25.51320775 -35.73124575   2.27482359 -23.67720440   7.67981459 -30.63388731
 [883]  -2.12532769  -6.06248123 -14.67967251   2.92695069 -32.55242308 -19.80182640 -12.70340060
 [890]  -0.36473422 -24.33804299 -33.53505272   6.38777505  -7.65940679 -45.22813407  -5.91512961
 [897]   4.82722129 -32.55034135  -5.64002137  -1.85377692 -24.82298659 -11.91899896   5.90226019
 [904]  -6.67799556 -18.39929702 -14.70709248  20.16465530  -2.37785503 -19.40544013 -43.64259489
 [911]  -4.80310727 -20.67267587  -6.12960286  14.76051916 -24.94995895 -21.55367734   3.51347606
 [918] -21.82098554   4.68892318  22.32281743   3.01554647 -14.22391287 -28.44488042   0.32000549
 [925] -26.29548705 -28.39677088  -6.06084948  -2.71491964   1.69227810  -8.71016310  -6.16547536
 [932] -11.67413566 -11.59680714   9.40984647  -9.93669428 -17.84745893   2.04601218 -19.80104095
 [939]  -0.49341925   6.14760676 -22.21183010  13.50485022  -5.22057307 -17.82539558 -24.46518962
 [946]  13.48666595 -11.80484500 -20.85173165 -31.08852302  -4.75860232 -36.88054918   2.98370161
 [953]  -0.37405972  17.58168372   0.48972051   9.79846880 -45.38152568 -25.01098967  -0.10839639
 [960]   0.14270290  10.28628008  21.54101027 -17.49673999 -12.24470665 -11.99999059 -24.05651849
 [967] -18.49661342 -12.30226946  -6.43942483 -23.07187196  -1.29013458 -15.53084359   0.86159642
 [974] -11.26368153  -4.91450784  -9.73711163 -14.70304307  -4.39913543 -21.76712550  -7.49898974
 [981] -25.17421015 -10.35273557  -9.64669400 -14.19622872 -13.91668603 -24.13258717 -15.08519499
 [988]   1.35746984  10.40157841  -2.47480562 -35.30199191 -25.52554695  -2.31850569 -10.24616931
 [995] -22.27223290  -4.57167529  -7.75456863  -2.13869306 -17.30982789 -24.04030147
> t.cal.rm <- mean(diff.s)/se(diff.s)
> t.cal.rm
[1] -19.82846
> p.val.rm <- pt(abs(t.cal.rm), length(s1)-1, lower.tail = F) * 2
> p.val.rm
[1] 4.836389e-74
> t.test(s1, s2, paired = T)

	Paired t-test

data:  s1 and s2
t = -19.828, df = 999, p-value < 2.2e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -10.073732  -8.259379
sample estimates:
mean difference 
      -9.166555 

> 
> # create multiple histogram
> s.all <- c(s1,s2)
> mean(s.all)
[1] 105.2971
> hist(s1, col='grey', breaks=50, xlim=c(50, 150))
> hist(s2, col='darkgreen', breaks=50, add=TRUE)
> abline(v=c(mean(s.all)), 
+        col=c("red"), lwd=3)
> abline(v=c(mean(s1), mean(s2)), 
+        col=c("black", "green"), lwd=3)
> 
> comb <- data.frame(s1,s2)
> dat <- stack(comb)
> head(dat)
     values ind
1  93.17788  s1
2 103.00254  s1
3 104.53388  s1
4  88.59698  s1
5 105.67789  s1
6 112.72657  s1
> 
> m.tot <- mean(s.all)
> m.s1 <- mean(s1)
> m.s2 <- mean(s2)
> 
> ss.tot <- ss(s.all)
> ss.s1 <- ss(s1)
> ss.s2 <- ss(s2)
> 
> df.tot <- length(s.all)-1
> df.s1 <- length(s1)-1
> df.s2 <- length(s2)-1
> 
> ms.tot <- var(s.all)
> ms.tot
[1] 125.5892
> ss.tot/df.tot
[1] 125.5892
> 
> var(s1)
[1] 108.3973
> ss.s1 / df.s1
[1] 108.3973
> 
> var(s2)
[1] 100.8519
> ss.s2 / df.s2
[1] 100.8519
> 
> ss.b.s1 <- length(s1) * ((m.tot - m.s1)^2)
> ss.b.s2 <- length(s2) * ((m.tot - m.s1)^2)
> ss.bet <- ss.b.s1+ss.b.s2
> ss.bet
[1] 42012.87
> 
> ss.wit <- ss.s1 + ss.s2
> ss.wit
[1] 209039.9
> 
> ss.bet + ss.wit
[1] 251052.8
> ss.tot
[1] 251052.8
> 
> library(dplyr)

다음의 패키지를 부착합니다: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union
> # df.bet <- length(unique(dat)) - 1
> df.bet <- nlevels(dat$ind) - 1
> df.wit <- df.s1+df.s2
> df.bet
[1] 1
> df.wit
[1] 1998
> df.bet+df.wit
[1] 1999
> df.tot 
[1] 1999
> 
> ms.bet <- ss.bet / df.bet
> ms.wit <- ss.wit / df.wit
> ms.bet
[1] 42012.87
> ms.wit
[1] 104.6246
> 
> f.cal <- ms.bet / ms.wit
> f.cal
[1] 401.5582
> pf(f.cal, df.bet, df.wit, lower.tail = F)
[1] 1.522061e-81
> 
> 
> f.test <-  aov(dat$values~ dat$ind, data = dat)
> summary(f.test)
              Df Sum Sq Mean Sq F value Pr(>F)    
dat$ind        1  42013   42013   401.6 <2e-16 ***
Residuals   1998 209040     105                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> sqrt(f.cal)
[1] 20.03892
> t.cal.ts
[1] -20.03892
> 
> # this is anova after all. 
> 
</code>