# 迴圈練習：*T* 檢定與中央及限定裡模擬
###### tags: `R教學` `基礎教學`

rm(list=ls())

#### t test p-value calculation ####
t.test(x, mu = 7)
t.test(rnorm(20, 1, 1), rnorm (20, 2, 1))

#### normal distribution and mean ####
curve(dnorm(x, 5, 2),lwd = 3, xlim = c(-5,15))
x <- rnorm(10, mean = 5, sd = 2)
mean(x)

##### 中央極限定理 #####
##1. 從常態分佈 (mean = 5, sd = 2) 中隨機抽取 10 個樣本，其平均的分布
mean (rnorm(10, mean = 5, sd = 2))

c.m <- numeric(1000)
c.sd <- numeric(1000)
for( i in 1:1000){
  a <- rnorm(10, 5, 2)
  c.m[i] <- mean(a)
  c.sd[i] <- sd(a)
}

hist(c.m, freq=F, breaks = 19, main = "樣本數 = 10")
curve(dnorm(x, 5, 2/sqrt(10)), add=T, col = 1, lty = 1, lwd = 2)

##2. 把平均的分布標準化
s.c.m <- (c.m-5)/(c.sd/sqrt(10))
hist(s.c.m, freq=F, ylim= c(0,0.4),  main = "樣本數 = 10", breaks = 19)
curve(dnorm(x, 0, 1), add=T, col = 1, lty = 1, lwd = 2)

polygon (c(seq(-4, qnorm(0.025, 0, 1), length = 100), rev(seq(-4, qnorm(0.025, 0, 1), length = 100))),
         c(dnorm (seq(-4, qnorm(0.025, 0, 1), length = 100), 0, 1),rep(0, 100)),
         col = rgb(red=1, green=0, blue=0, alpha = 0.5), border = NA)
polygon (c(seq(qnorm(0.975, 0, 1), 4, length = 100), rev(seq(qnorm(0.975, 0, 1), 4, length = 100))),
         c(dnorm (seq(qnorm(0.975, 0, 1), 4, length = 100), 0, 1),rep(0, 100)),
         col = rgb(red=1, green=0, blue=0, alpha = 0.5), border = NA)

##3. 信賴區間 
text(-2, 0.15, paste0("down = ",round(qnorm(0.025, 0, 1)*(2/sqrt(10))+5, 2)), cex = 1.5)
text( 2, 0.15, paste0("up = ", round(qnorm(0.975, 0, 1)*(2/sqrt(10))+5, 2)), cex = 1.5)


##4. 理論數值與實際數值
length(s.c.m[s.c.m<qnorm(0.025, 0, 1)|
             s.c.m>qnorm(0.975, 0, 1)])/1000

##5. T 分布 and 常態分佈
hist(s.c.m, freq=F, ylim= c(0,0.4),  main = "樣本數 = 10", breaks = 19)
curve(dnorm(x, 0, 1), add=T, col = 1, lty = 1, lwd = 2)
curve(dt(x, 9), add=T, col = 2, lty = 1, lwd = 2)

abline (v = qnorm(0.025, 0, 1), lwd = 2)
abline (v = qt(0.025,9), col = 2, lwd = 2)
#根據常態切位的實際數值
text(-1.8, 0.35, paste0("PN = ", length(s.c.m[s.c.m<qnorm(0.025, 0, 1)|s.c.m>qnorm(0.975, 0, 1)])/1000), 
     adj=c(0,1), cex = 1.5)
text(-3.5, 0.35, paste0("PT = ", length(s.c.m[s.c.m<qt(0.025, 9)|s.c.m>qt(0.975, 9)])/1000), 
     adj=c(0,1), cex = 1.5, col = 2)

length(s.c.m[s.c.m<qnorm(0.025, 0, 1)|
               s.c.m>qnorm(0.975, 0, 1)])/1000
length(s.c.m[s.c.m<qt(0.025, 9)|
               s.c.m>qt(0.975, 9)])/1000

t.test(x, mu = 7)
(mean(x)-7)/(sd(x)/sqrt(10))
pt(-4.976963, 9)*2