Data Loading

load("./crimtab.RData")
ls()
## [1] "crimtab.2"       "crimtab.df"      "crimtab.long"    "crimtab.long.df"
ls.str()
## crimtab.2 :  'table' int [1:42, 1:22] 0 0 0 0 0 0 1 0 0 0 ...
## crimtab.df : 'data.frame':   924 obs. of  3 variables:
##  $ finger: num  9.4 9.5 9.6 9.7 9.8 9.9 10 10.1 10.2 10.3 ...
##  $ height: num  56 56 56 56 56 56 56 56 56 56 ...
##  $ Freq  : int  0 0 0 0 0 0 1 0 0 0 ...
## crimtab.long :  num [1:3000, 1:2] 10 10.3 9.9 10.2 10.2 10.3 10.4 10.7 10 10.1 ...
## crimtab.long.df : 'data.frame':  3000 obs. of  2 variables:
##  $ finger: num  10 10.3 9.9 10.2 10.2 10.3 10.4 10.7 10 10.1 ...
##  $ height: num  56 57 58 58 58 58 58 58 59 59 ...
head(crimtab.long.df, n = 10)
##    finger height
## 1    10.0     56
## 2    10.3     57
## 3     9.9     58
## 4    10.2     58
## 5    10.2     58
## 6    10.3     58
## 7    10.4     58
## 8    10.7     58
## 9    10.0     59
## 10   10.1     59

Student 의 Simulation 재현

3,000장의 카드를 잘 섞는 것은 sample() 이용.

# set.seed(113)
crimtab.shuffle <- crimtab.long.df[sample(1:3000), ]
head(crimtab.shuffle, n = 10)
##      finger height
## 1273   11.5     65
## 2192   11.8     67
## 1580   11.2     66
## 498    11.1     63
## 1439   11.9     65
## 420    10.9     63
## 2463   11.6     68
## 2190   11.8     67
## 806    11.2     64
## 471    11.0     63

표본의 크기가 4인 750개의 표본을 만드는 작업은 rep() 이용.

sample.id <- as.factor(rep(1:750, each = 4))
head(sample.id, n = 10)
##  [1] 1 1 1 1 2 2 2 2 3 3
## 750 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 ... 750

각 표본의 평균과 표준편차 계산에는 tapply() 이용.

finger.sample.mean <- tapply(crimtab.shuffle[, "finger"], sample.id, mean)
finger.sample.sd <- tapply(crimtab.shuffle[, "finger"], sample.id, sd)
str(finger.sample.mean)
##  num [1:750(1d)] 11.4 11.6 11.2 11.9 12 ...
##  - attr(*, "dimnames")=List of 1
##   ..$ : chr [1:750] "1" "2" "3" "4" ...
str(finger.sample.sd)
##  num [1:750(1d)] 0.316 0.451 0.126 0.591 1.164 ...
##  - attr(*, "dimnames")=List of 1
##   ..$ : chr [1:750] "1" "2" "3" "4" ...

t-통계량 계산. Student는 표준편차 계산에서 분모에 \(n\)을 사용하고 히스토그램을 그려 비교하였으나 자유도 3인 t-분포와 비교하기 위하여 \(t=\frac{\bar{X_n}-\mu}{\hat{SD}/\sqrt{n}}\)을 계산함. (여기서 \(\hat{SD}\)는 표본 표준편차)

sample.t <- (finger.sample.mean - mean(crimtab.long.df[, "finger"]))/(finger.sample.sd/sqrt(4))
str(sample.t)
##  num [1:750(1d)] -0.932 0.0117 -5.9185 1.2782 0.8204 ...
##  - attr(*, "dimnames")=List of 1
##   ..$ : chr [1:750] "1" "2" "3" "4" ...

계산한 t-통계량 값들의 평균과 표준편차, 히스토그램을 그리고 자유도 3인 t-분포의 밀도함수 및 표준정규곡선과 비교. 우선 모두 같은 값들이 나와서 분모가 0인 경우가 있는지 파악. 있으면 모평균과 비교하여 양수인 경우 +6, 음수인 경우 -6 값 부여(Student가 한 일)

t.inf <- is.infinite(sample.t)
sample.t[t.inf]
## named numeric(0)
sample.t[t.inf] <- 6*sign(sample.t[t.inf])

문제되는 값이 없는 것을 확인하고, 평균과 표준편차 계산. 자유도 \(n\)인 t-분포의 평균과 표준편차는 각각 0과 \(\sqrt{\frac{n}{n-2}}\)임을 상기할 것. -6이나 +6보다 큰 값이 상당히 자주 나온다는 점에 유의.

mean(sample.t)
## [1] -0.02539319
sd(sample.t)
## [1] 1.558049
summary(sample.t)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -7.41800 -0.80190  0.01053 -0.02539  0.86230  7.14900

t-통계량들의 히스토그램을 그리고, 자유도 3인 t의 밀도함수, 표준정규분포 밀도함수와 비교.

# hist(sample.t, prob = TRUE, ylim = c(0, 0.5))
# hist(sample.t, prob = TRUE, nclass = 20, xlim = c(-6, 6), ylim = c(0, 0.5), main = "Histogram of Sample t-statistics", xlab = "Sampled t-values")
# hist(sample.t, prob = TRUE, nclass = 50, xlim = c(-6, 6), ylim = c(0, 0.5), main = "Histogram of Sample t-statistics", xlab = "Sampled t-values")
hist(sample.t, prob = TRUE, breaks = seq(-20, 20, by = 0.5), xlim = c(-6, 6), ylim = c(0, 0.5), main = "Histogram of Sample t-statistics", xlab = "Sampled t-values")
lines(seq(-6, 6, by = 0.01), dt(seq(-6, 6, by = 0.01), df = 3), col = "blue")
lines(seq(-6, 6, by = 0.01), dnorm(seq(-6, 6, by = 0.01)), col = "red")
legend("topright", inset = 0.05, lty = 1, col = c("blue", "red"), legend = c("t with df = 3", "standard normal"))

qqnorm() 을 그려보면 정규분포와 꼬리에서 큰 차이가 난다는 것을 알 수 있음.

qqnorm(sample.t)