> 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] 194
> 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, 49))

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

tout <- t.test(height_sample_biased, mu=pop_mean)
tout

> 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 <- 36
s.mean <- 53.5
s.sd <- 5
set.seed(1)
s1 <- rnorm(s.n, s.mean, s.sd) 
s1
summary(s1)

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

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

t.test(s1, mu=51)

표준편차가 샘플과 모집단이 서로 다르다고 한다면

t.test(s1, mu=51, var.equal=FALSE)

# 이것을 손으로 풀면 
# 가설: 강사의 학습방법이 효과가 있을 것이다
# 가설: 강사의 학생들의 점수와 모집단의 평균점수가 다를 것이다. 

# 모집단 평균: 51
s.mu <- 51
# 강사집단의 평균: 53.5, sd = 5

# 이 때 se 값은 
se <- sqrt(s.sd^2/s.n)
se2 <- se * 2
# 영가설 기각 범위는 아래 두 점수 사이
s.mu - se2
s.mu + se2
# mean(s1)값이 이 점수 밖에 존재하는가?
mean(s1)

t-test

게임회사의 직원인 강태원씨는 게임에 접속하는 사용자의 게임에 대한 강박성이 존재할 때, 가장 접속을 많이 할 것이라는 가설을 세웠다. 강박성이란 게임의 진행과정이 머리에 끊임없이 떠올라서 게임상황을 활성화하는 정도를 말한다. 이런 방법으로 강태원은 강박성이 높은 그룹과 낮은 그룹 두 그룹을 구하고, 각 그룹의 게임에 쓴 비용을 구하였다. 정말 차이가 있을까?

HI group
평균: 6000
n: 16
LO group
평균: 4000
n: 16

s.n <- 16
s.hi <- 6000
s.sd <- 500
s.lo <- 4000

set.seed(101)
s.hi <- round(rnorm(s.n, s.hi, s.sd))
s.lo <- round(rnorm(s.n, s.lo, s.sd))

s.hi
s.lo

mean(s.hi)
mean(s.lo)

t.test(s.hi, s.lo, var.equal=T)

미디어교육을 초등학교 5학년 학생에게 하여, 이 학생들이 게임에 쓴 시간을 조절을 할 수 있게 된다고 학부모들을 설득하려고 한다. 이를 위해서 강경태 게임회사 기획자는 미디어교육 방법을 개발하고 이를 5학년 학생 집단 10명에게 적용한 후 게임에 쓴 시간을 측정을 하였다. 그 결과는 아래와 같다.

• se 값을 구해보시오
• t 값을 구해보시오

before <- c(392.9, 393.2, 345.1, 393, 434, 427.9, 422, 383.9, 392.3, 352.2)
after <- c(200.1, 190.9, 192.7, 213, 241.4, 196.9, 172.2, 185.5, 205.2, 193.7)

bminusa <- before-after
mean(bminusa)
sd(bminusa)

se <- sqrt(sd(bminusa)^2/10)
t.value <- (mean(before)-mean(after))/se

t.test(before, after, paired=TRUE)