> rnorm2 <- function(n,mean,sd) { mean+sd*scale(rnorm(n)) }
> potato_sample <- rnorm2(25, 194,20)
> mean(potato_sample)
[1] 194
> sqrt(var(potato_sample))
     [,1]
[1,]   20

> t.test(potato_sample, mu=200)

	One Sample t-test

data:  potato_sample
t = -1.5, df = 24, p-value = 0.1467
alternative hypothesis: true mean is not equal to 200
95 percent confidence interval:
 185.7444 202.2556
sample estimates:
mean of x 
      194 
>       
> abs(qt(0.05/2, 24))
[1] 2.063899

> rnorm2 <- function(n,mean,sd) { mean+sd*scale(rnorm(n)) }
> potato_sample_large <- rnorm2(2500, 194,20)
> mean(potato_sample_large)
[1] 191
> sqrt(var(potato_sample_large))
[1]   20
> t.test(potato_sample_large, mu=200)

	One Sample t-test

data:  potato_sample_large
t = -15, df = 2499, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 200
95 percent confidence interval:
 193.2156 194.7844
sample estimates:
mean of x 
      194 

> abs(qt(0.05/2, 2499))
[1] 1.960914      
> 

Also refer to http://stats.seandolinar.com/one-sample-t-test-with-r-code/

#reads data set
data2 <- read.csv("http://commres.net/wiki/_media/r/Height_data2.csv", header=T)
height <- data2$height

#N - number in population
#n - number in sample
N <- length(height)
n <- 50

#population mean
pop_mean <- mean(height)
pop_sd <- sd(height)

#tall-biased sample
cut <- 1:25000
weights <- cut^.6
sorted_height <- sort(height)
set.seed(123)
height_sample_biased <- sample(sorted_height, size=n)

sample_mean <- mean(height_sample_biased)
sample_sd <- sd(height_sample_biased)

#t-stat
t <- (sample_mean - pop_mean) / (sample_sd/sqrt(n))

abs(qt(0.05/2, 50))

# or 
1-pt(t,n-1)
# n = 50 이었음
# t = 계산된 값

> rnorm2 <- function(n,mean,sd) { mean+sd*scale(rnorm(n)) }
> bird_eye <- rnorm2(16, 35, 9)
> hist(bird_eye)
> bird_eye
          [,1]
 [1,] 33.42614
 [2,] 34.29464
 [3,] 40.73375
 [4,] 38.16267
 [5,] 39.88801
 [6,] 46.02066
 [7,] 48.00472
 [8,] 41.13903
 [9,] 20.64618
[10,] 32.34706
[11,] 21.25446
[12,] 27.80643
[13,] 26.59595
[14,] 36.57292
[15,] 48.49734
[16,] 24.61004
attr(,"scaled:center")
[1] -0.1142795
attr(,"scaled:scale")
[1] 1.251198

> t.test(bird_eye, mu=30)

	One Sample t-test

data:  bird_eye
t = 2.2222, df = 15, p-value = 0.04207
alternative hypothesis: true mean is not equal to 30
95 percent confidence interval:
 30.20424 39.79576
sample estimates:
mean of x 
       35 

> abs(qt(0.05/2, 15))
[1] 2.13145

> rnorm2 <- function(n,mean,sd) { mean+sd*scale(rnorm(n)) }
> groupA <- rnorm2(10, 19, sqrt(160/9))
> groupB <- rnorm2(10, 25, sqrt(200/9))
> t.test(groupA, groupB)

	Welch Two Sample t-test

data:  groupA and groupB
t = -3, df = 17.78, p-value = 0.007762
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10.205566  -1.794434
sample estimates:
mean of x mean of y 
       19        25 

> t.test(groupA, groupB, var.equal=T)

	Two Sample t-test

data:  groupA and groupB
t = -3, df = 18, p-value = 0.007685
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10.201844  -1.798156
sample estimates:
mean of x mean of y 
       19        25 

> abs(qt(.05/2, 18))
[1] 2.100922
> 

> data(anorexia, package="MASS")       # weight gain (lbs.) in anorexic women
> attach(anorexia)
> str(anorexia)
'data.frame':	72 obs. of  3 variables:
 $ Treat : Factor w/ 3 levels "CBT","Cont","FT": 2 2 2 2 2 2 2 2 2 2 ...
 $ Prewt : num  80.7 89.4 91.8 74 78.1 88.3 87.3 75.1 80.6 78.4 ...
 $ Postwt: num  80.2 80.1 86.4 86.3 76.1 78.1 75.1 86.7 73.5 84.6 ...
> ft = subset(anorexia, subset=(Treat=="FT"))         # just the family therapy threatment
> ft
   Treat Prewt Postwt
56    FT  83.8   95.2
57    FT  83.3   94.3
58    FT  86.0   91.5
59    FT  82.5   91.9
60    FT  86.7  100.3
61    FT  79.6   76.7
62    FT  76.9   76.8
63    FT  94.2  101.6
64    FT  73.4   94.9
65    FT  80.5   75.2
66    FT  81.6   77.8
67    FT  82.1   95.5
68    FT  77.6   90.7
69    FT  83.5   92.5
70    FT  89.9   93.8
71    FT  86.0   91.7
72    FT  87.3   98.0
> detach(anorexia)
> rm(anorexia)
> t.test(ft$Prewt, ft$Postwt, paired=T)

	Paired t-test

data:  ft$Prewt and ft$Postwt
t = -4.1849, df = 16, p-value = 0.0007003
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10.94471  -3.58470
sample estimates:
mean of the differences 
              -7.264706 

> 

> t.test(ft$Postwt-ft$Prewt, mu=0, data=ft)

	One Sample t-test

data:  ft$Postwt - ft$Prewt
t = 4.1849, df = 16, p-value = 0.0007003
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
  3.58470 10.94471
sample estimates:
mean of x 
 7.264706 

>

내가 키운 감자하나의 평균은 50g 이고 표준편차는 4g 이다. 이 감자에서 n = 25 샘플을 취하여 평균을 구하고 이를 반복하여 1000번 해보자.

• 이 때, 샘플의 평균은 어떤 분포를 이루어야 하는가?
• 이 집합의 표준편차값은 무엇인가?
• 강사가 (CLT 설명에서) 제시한 계산값과 비교해보라.
• 표준편차 2개를 써서 평균에 더하고 뺀 값의 밖에 존재하는 갯수는 몇개인지 세어보라

n <- 25
s.mean <- 50
s.sd <- 4
s1 <- rnorm(s.size, s.mean, s.sd) 
s1
summary(s1)

그런데, 문제는 이를 100번 반복해 보고 이 평균을 기록하여 기록된 것의 분포를 보라는 이야기이다.

iter <- 1000
mean.dist <- rep (NA, iter)
s.n <- 100
s.mean <- 50
s.sd <- 4
for (i in 1:iter)  {
  mean.dist[i] <- mean(rnorm(s.n, s.mean, s.sd))
}
summary(mean.dist)
mean(mean.dist) 
se.obt <- sd(mean.dist)
se <- sqrt((s.sd^2)/s.n)
se.obt
se # sd(mean.dist)값과 거의 같아야 함

se2 <- 2*se
s.mean 
se2

s.mean - se2
s.mean + se2

mean.dist > s.mean + se2
table(mean.dist > s.mean + se2)["TRUE"]

(mean.dist > s.mean + se2) | (mean.dist < s.mean - se2)
table((mean.dist > s.mean + se2) | (mean.dist < s.mean - se2))["TRUE"]

s.n <- 25
s.mean <- 52
s.sd <- 5
set.seed(1)
s1 <- rnorm(s.n, s.mean, s.sd) 
s1
summary(s1)

위의 샘플은 수학 강사인 A가 특별한 방법으로 학생들을 가르친 후 얻은 평균 점수이다 (52점, n=25). 전국평균이 50점임을 알고 있다고 할 때, 강사는 이 55점으로 자신의 학습방법이 효과가 있다는 것을 증명하려고 한다.

• 모집단의 표준오차값은 무엇인가?
• 자신의 학습방법이 효과가 있다는 것을 증명하시오
  • 가설
  • 영가설
  • 표준오차
  • 영가설 기각범위
  • 확인

t.test(s1, mu=50)

<