1. Fifty four university students (male and female) agreed to participate in an experiment on problem solving. They were tested on the vocabulary. Does the scores follow a normal distribution? If yes, what is the slope and the intercept? What are estimates of the mean and variance of normal distribution?

54명의 대학생들이 단어 시험에 참여해 각각 점수를 얻었다. 이 점수의 분포가 정규분포를 따르는지를 확인하기 위해서 두 가지 접근을 취할 것이다. 하나는 Normal Probability Plot을 통해 눈으로 확인하는 것이고, 나머지 하나는 비모수적 방법을 통해 정규분포를 따르는지에 대한 가설을 검정하는 것이다.

Normal Probability Plot

먼저 Normal Probability Plot을 그려보자. R은 이를 위한 qqnorm()이라는 함수를 제공한다. 결과는 아래와 같다.

Normal Probablity Plot의 점들이 직선을 이룬다면 우리는 이 분포가 정규분포를 따른다고 판단할 수 있다. 위의 그림은 눈으로 보기에 직선을 잘 이루고 있는 것처럼 보인다. 하지만 눈으로 판단하는 것에는 한계가 있으니 위 그림에 보조 직선을 그려보자. 보조 직선을 그리는 방법에는 두 가지 방법이 있다. 하나는 qqline()함수를 이용하는 것이고, 나머지 하나는 EDA적 방법을 통해 직선을 그리는 것이다. 두 가지 방법을 설명하기 위해 아래의 그림에는 plot을 생성하기 위한 코드를 그대로 첨부하였다.
먼저, qqline()이다. qqline()이라는 함수를 통해 간단하게 직선을 그릴 수 있지만 그 과정을 이해하기 위해 직선을 그리는 과정을 코드로 하나하나 작성해보았다. 제1사분위수와 제3사분위수를 관통하는 직선을 그리는 방법이다. 아래의 코드와 같은 방법으로 기울기와 절편을 각각 slope와 int라는 변수에 저장하였다. 그리고 초록색 선으로 직선을 그렸다.
두 번째 방법은 EDA적 방법으로 직선을 그렸다. 절편으로는 자료의 중위값을 사용하였으며, pseudosigma를 계산하여 이를 기울기로 사용하였다. EDA 직선은 빨간색으로 작성하였다.

qqnorm(score, ylab="Scores of Students",main="Normal Prob Plot")
#qqline(score)

score.x = quantile(score, probs=c(0.25,0.75), names = FALSE,na.rm = TRUE)
score.y = qnorm(c(0.25,0.75))
  
slope = diff(score.x)/diff(score.y)
int = score.x[1] - slope * score.y[1]
  
abline(int, slope, col="green2", lwd=3)

x = fivenum(score)
pseudosigma = (x[4]-x[2])/1.34
abline(x[3],pseudosigma, col="red", lwd=3)

legend(-2, 18, c("qqline", "EDA line"), col=c("green2", "red"), lty=1, lwd=3)


두 직선을 비교해 보았을 때, 빨간색 직선이 조금 더 자료를 잘 나타낸다고 볼 수 있다. 가장 위와 아래에 있는, 직선에서 약간 떨어진 두 점을 고려하면 EDA적으로 그린 빨간색 직선이 더 적합하다. 두 직선의 기울기는 거의 같아보이며 빨간색 직선이 약간 더 위로 수직이동한 모습을 확인할 수 있다. 실제로 두 직선의 기울기는 매우 비슷하며 절편에만 약간의 차이가 발생함을 수치로 확인할 수 있다. 초록선의 기울기는 2.223903이며 빨간선의 기울기는 2.238806이다. 초록선의 절편은 12.5, 빨간선의 절편은 13으로 0.5의 차이가 발생했다.

slope #초록선의 기울기
## [1] 2.223903
pseudosigma #빨간선의 기울기
## [1] 2.238806
int #초록선의 절편
## [1] 12.5
x[3] #빨간선의 절편
## [1] 13

여기서는 빨간선을 선택했으므로 Normal Prob Plot의 직선은 y=2.238806x + 13이라고 할 수 있다. 결론적으로, 54명의 대학생들이 단어 시험 점수 분포는 평균이 13이고 표준편차가 2.238806인 정규분포를 따른다고 말할 수 있다.

추가적으로 위 자료의 왜도를 살펴보았다. 왜도는 Hinge를 기준으로 EDA적으로 계산할 수 있으며 -1과 1의 범위를 가지고 0에 가까울수록 대칭이라고 할 수 있다. 대칭이라는 기준이 반드시 정규분포를 따른다고는 할 수 없지만 분포의 대칭성을 파악하는 일은 중요하다. 아래와 같이 왜도를 계산한 결과, -0.333이 도출되었다. 따라서 위의 자료는 약간의 대칭성을 가진다고 볼 수 있다. 물론 이 결과가 해당 자료의 정규분포 모양을 보장하는 것은 절대 아니다.

#skewness
skewness = function(x){
  Hl <- fivenum(x)[2]
  M <- fivenum(x)[3]
  Hu <- fivenum(x)[4]
  return((Hu-M-M+Hl)/(Hu-M+M-Hl))
}

skewness(score)
## [1] -0.3333333


또한, density graph를 통해서도 두 자료의 분포를 비교할 수 있다. 아래의 그림에서 검은색 선은 자료의 density distribution를 나타낸다. 빨간 선은 정규분포의 density distribution을 나타낸다. 단 여기서 빨간선을 그리는 정규분포는 위에서 우리가 선택한 평균이 13이고 표준편차가 2.238806인 정규분포이다. 한눈에 보기에도 자료의 분포가 정규분포와 굉장히 유사함을 확인할 수 있다.

비모수적 검정

분포의 모양을 비교하는 비모수적 검정에는 크게 3가지가 있다.

  • Kolmogorov-Smirnov test
  • Cramer-von Mises test
  • Anderson-Darling test

자료의 분포 비교를 해주는 검정인데 공통적으로 귀무가설은 “두 자료의 CDF가 같다.”이다. 일반적으로 Kolmogorov-Smirnov test를 자주 사용하지만 검정의 power 측면에서는 Anderson-Darling test가 가장 강력한 것으로 알려져 있다. 위의 검정들을 통해 54명의 대학생들이 단어 시험 점수 분포가 정규분포를 따르는지를 검정해보자. 여기서 기준이 되는 정규분포는 우리가 위에서 선택한 평균이 13이고 표준편차가 2.238806인 정규분포이다.

#Kolmogorov-Smirnov test
ks.test(score, "pnorm", x[3], pseudosigma)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  score
## D = 0.14749, p-value = 0.1907
## alternative hypothesis: two-sided
#Cramer-von Mises test
library(goftest)
cvm.test(score, "pnorm", x[3], pseudosigma)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  Parameters assumed to be fixed
## 
## data:  score
## omega2 = 0.20745, p-value = 0.2539
#Anderson-Darling test
ad.test(score, "pnorm", x[3], pseudosigma)
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  Parameters assumed to be fixed
## 
## data:  score
## An = 1.1983, p-value = 0.268

KS검정의 p값은 0.1907, CVM검정의 p값은 0.2539, AD검정의 p값은 0.268이다. 유의수준 0.05를 기준으로 할 때 모든 검정에서 귀무가설을 기각할 수 없다. 따라서, 54명의 대학생들이 단어 시험 점수 분포는 평균이 13이고 표준편차가 2.238806인 정규분포를 따른다고 말할 수 있다.







2. In 1960s in the United States of America 16 states owned the retail liquor stores while in 26 the stores were privately owned. (Some are omitted for technical reasons.) The table shows the price in dollars of a fifth of Seagram 7 Crown Whisky in two sets of states in 1961. Do the distributions of prices of two group follow the same distribution?


QQ-plot

두 자료의 비교 역시 qqplot()을 통해 그림을 그리고 값들이 직선 상에 있는지를 통해 확인할 수 있다. 다만 두 자료의 크기가 다른 경우에는 크기가 큰 자료의 값들이 작은 자료의 크기에 맞추어 interpolated되는 것을 유념하자. 실제로 monopoly와 private에 대해 qqplot을 그리면 값들은 16개가 그려지는 것을 확인할 수 있다.

monopoly = c(4.65,  4.55,  4.11,  4.15,  4.20,  4.55,  3.80,  4.00, 4.19,  4.75,
             4.74,  4.50,  4.10,  4.00,  5.05,  4.20)
private = c(4.82,  5.29,  4.89,  4.95,  4.55,  4.90,  5.25,  5.30,  4.29,  4.85,
            4.54,  4.75,  4.85,  4.85,  4.50,  4.75,  4.79,  4.85,  4.79,  4.95,
            4.95,  4.75,  5.20,  5.10,  4.80,  4.29)
length(monopoly)
## [1] 16
length(private)
## [1] 26

qqplot()의 결과는 아래와 같다. 그림에 qqline()도 같이 그려주었다. 전체적으로 값들이 직선을 이루고 있기 때문에 두 자료의 분포는 같은 모양이라고 말할 수 있다. 즉 monopoly와 private한 상점에서의 위스키 가격의 분포는 같은 모양이라고 말할 수 있다.

## 
## Call:
## line(qqplot(monopoly, private))
## 
## Coefficients:
## [1]  2.1013  0.6382


Tukey Mean-Difference Plot

위에서 두 자료의 분포가 동일한 모양을 갖는다고 했지만 실제로 값들을 살펴보면 약간의 차이가 있음을 확인할 수 있다. 가령 두 분포가 동일한 모양을 가진다고 하더라도 위치에 차이가 있을 수 있다. 이를 확인하기 위해 Tukey의 Mean-Difference Plot을 그려보았다.
점들이 수평으로 직선을 이룰 때가 가장 이상적인 모습이다. 다만 위 그림에서 자료의 값들은 수평으로 직선을 이루지만 0에서가 아니라 0.5에서 직선을 이룬다. 즉 Monopoly와 Private에 0.5의 차이가 평균적으로 발생한다는 의미이다. 분포의 모양은 동일하지만 평균적으로 Private 상점의 위스키 가격이 약 0.5달러 정도 더 비싸다는 것이다.

Kolmogorov-Smirnov test

위의 결과를 가설검정하기 위해 KS검정을 진행했다. 귀무가설은 “monopoly와 private의 위스키 가격의 CDF는 같다.”이며 양측검정이다. 결과적으로 p값이 0.0001956이 나와 유의수준 0.05를 기준으로 귀무가설을 기각했다. 즉 두 분포의 CDF는 다르다.

ks.test(monopoly, private, alternative="two.sided")
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  monopoly and private
## D = 0.68269, p-value = 0.0001956
## alternative hypothesis: two-sided

이번에는 Tukey’s Mean-Difference Plot의 결과해석을 반영하여 monoploy에 0.5를 더하여 KS검정을 해보았다. p값은 0.1324로 귀무가설을 기각할 수 없다. 즉 두 분포의 CDF는 같다고 말할 수 있다.

ks.test(monopoly+0.5, private, alternative="two.sided")
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  monopoly + 0.5 and private
## D = 0.37019, p-value = 0.1324
## alternative hypothesis: two-sided

위의 모든 결과를 종합하면 monopoly와 private 상점에서의 위스키 가격의 분포 모양은 동일하지만 private의 자료가 monopoly보다 0.5정도 오른쪽으로 이동해있다고 해석할 수 있다.

ROC curve and AUC

추가적으로 ROC curve와 AUC를 그려보자.

ROC 곡선의 해석은 가운데 대각선으로부터 곡선이 얼마나 떨어져있는가로 판단한다. 곡선이 가운데 대각선으로부터 멀리 떨어져 있을수록 두 자료의 분포는 다르다고 말할 수 있다. 왼쪽의 그림은 스케일링을 하지 않은 그대로의 monopoly와 private의 ROC 곡선이다. 곡선이 대각선으로부터 멀리 떨어져있다. 오른쪽은 monopoly에 0.5를 더해주어 스케일링을 해준 것이다. 오른쪽의 ROC 곡선은 대각선을 따라서 그려졌음을 확인할 수 있다. 이를 수치적으로 확인할 수도 있다. AUC(Area Under Curve)라고 불리는 통계량은 ROC 곡선의 아래의 면적을 나타내는 수치이다. 이 값이 1에 가까울수록 두 자료의 분포는 다르다고 말할 수 있으며 0.5에 가까울수록 두 자료의 분포는 같다고 말할 수 있다. 두 그림에 대해 AUC를 계산해보면 왼쪽은 0.8822, 오른쪽은 0.5325이다. 따라서 monopoly와 private 상점에서의 위스키 가격의 분포 모양은 동일하지만 private의 자료가 monopoly보다 0.5정도 오른쪽으로 이동해있다고 해석할 수 있다.

#왼쪽그림의 AUC
auc(res) 
## Area under the curve: 0.8822
#오른쪽그림의 AUC
auc(res.ad)
## Area under the curve: 0.5325















3. The frequency table gives the yield strength of a Bofores steel. The observed values are obtained as routine tests. x=yield strength in 1.275 kg/. Does the yield strength follow a normal distribution? Does the data fit to an exponential distribution? Does the data fit to Weibull distribution better? What are the estimates of the parameters of Weibull distribution?


누적도수분포표를 바탕으로 원자료를 도출해냈다. 원자료를 가지고 Normal Probability Plot, Exponential Probability Plot, Weibull Probability Plot을 그렸다. Exponential Prob Plot에서는 원점에 너무 많은 자료가 몰려있어서 1/3승하여 재표현을 해주는 그림을 추가로 그려주었다.

위의 4개의 그림을 종합적으로 살펴보면 Normal Prob Plot의 값들이 가장 직선을 이루고 있는 것을 확인할 수 있다. 기본적으로 qqplot은 값들이 직선을 이루어야 해당 분포를 따른다고 말할 수 있다. Normal을 제외한 Exponential과 Weibull의 Prob Plot은 곡선을 이루고 있기 때문에 위의 자료가 해당 분포들을 따른다고 말하기 어렵다. (세번째 그림인 exp에 1/3승을 해준 자료는 재표현을 해주었기 때문에 자료가 직선을 이룬다고 하더라도 Exponential 분포를 따른다고는 말할 수 없다.)

따라서 위 자료는 정규분포를 따른다고 말할 수 있다. 이에 대한 근거를 보충하기 위해 평균과 분산의 추정치를 구해 qqline()을 그리고 히스토그램을 그려보자. 먼저 qqline()을 그리기 위해 기울기와 절편을 계산한다. 기울기는 2.223903이며 절편은 36.5이다.

score.x = quantile(x, probs=c(0.25,0.75), names = FALSE,na.rm = TRUE)
score.y = qnorm(c(0.25,0.75))
  
slope = diff(score.x)/diff(score.y)
int = score.x[1] - slope * score.y[1]

slope
## [1] 2.223903
int
## [1] 36.5


그래프를 그린 결과 모든 점들이 직선을 잘 이루고 있음을 확인할 수 있다. 기울기는 2.223903, 절편은 36.5이므로 위의 자료는 평균 36.5, 표준편차 2.223903인 정규분포를 따른다고 할 수 있다.

Weibull 분포는 위의 그림을 살펴볼 때, 정규분포보다 적합해보이지는 않는다. 하지만 Weibull 분포의 parameter의 추정치를 구해보자. 이를 위해서는 qqline을 그려 기울기와 절편을 확인해야한다. 결과적으로 그린 직선의 기울기는 0.04876, 절편은 3.61148이다.
그림은 아래와 같다.
Weibull 분포의 parameter의 추정치는 아래의 수식으로 구할 수 있다.

(b=1/0.04876)
## [1] 20.50861
(a=exp(3.61148*(-b)))
## [1] 6.813188e-33

따라서, Weibull 분포의 shape parameter인 a와 scale parameter인 b는 각각 6.813188e-33과 20.50861로 추정할 수 있다.










R Script

보고서에 쓰인 모든 코드를 아래에 정리해두었다.

#1.
score = c(14, 11, 13, 13, 13, 15, 11, 16, 10,
          13, 14, 11, 13, 12, 10, 14, 10, 14,
          16, 14, 14, 11, 11, 11, 13, 12, 13,
          11, 11, 15, 14, 16, 12, 17,  9, 16,
          11, 19, 14, 12, 12, 10, 11, 12, 13,
          13, 14, 11, 11, 15, 12, 16, 15, 11)

qqnorm(score, ylab="Scores of Students",main="Normal Prob Plot")

qqnorm(score, ylab="Scores of Students",main="Normal Prob Plot")
#qqline(score)

score.x = quantile(score, probs=c(0.25,0.75), names = FALSE,na.rm = TRUE)
score.y = qnorm(c(0.25,0.75))

slope = diff(score.x)/diff(score.y)
int = score.x[1] - slope * score.y[1]

abline(int, slope, col="green2", lwd=3)

x = fivenum(score)
pseudosigma = (x[4]-x[2])/1.34
abline(x[3],pseudosigma, col="red", lwd=3)

legend(-2, 18, c("qqline", "EDA line"), col=c("green2", "red"), lty=1, lwd=3)

slope #초록선의 기울기
pseudosigma #빨간선의 기울기
int #초록선의 절편
x[3] #빨간선의 절편

#skewness
skewness = function(x){
  Hl <- fivenum(x)[2]
  M <- fivenum(x)[3]
  Hu <- fivenum(x)[4]
  return((Hu-M-M+Hl)/(Hu-M+M-Hl))
}

skewness(score)

#density function
hist(score, freq = FALSE, ylim=c(0,0.3))
lines(density(score), lwd=2)
n = rnorm(100000, x[3], pseudosigma)
lines(density(n), col="red", lwd=2)
legend(18, 0.28, c("Score", "Normal"), col=c("black", "red"), lty=1, lwd=2)

#Kolmogorov-Smirnov test
ks.test(score, "pnorm", x[3], pseudosigma)
#Cramer-von Mises test
library(goftest)
cvm.test(score, "pnorm", x[3], pseudosigma)
#Anderson-Darling test
ad.test(score, "pnorm", x[3], pseudosigma)


#2. 
monopoly = c(4.65,  4.55,  4.11,  4.15,  4.20,  4.55,  3.80,  4.00, 4.19,  4.75,
             4.74,  4.50,  4.10,  4.00,  5.05,  4.20)
private = c(4.82,  5.29,  4.89,  4.95,  4.55,  4.90,  5.25,  5.30,  4.29,  4.85,
            4.54,  4.75,  4.85,  4.85,  4.50,  4.75,  4.79,  4.85,  4.79,  4.95,
            4.95,  4.75,  5.20,  5.10,  4.80,  4.29)
length(monopoly)
length(private)

qqplot(monopoly, private)
line(qqplot(monopoly, private))
abline(line(qqplot(monopoly, private)))
title("qqplot ; monopoly states vs private-ownership states")


#tukey mean-diff
qq.x = qqplot(monopoly,private)$x
qq.y = qqplot(monopoly,private)$y
plot((qq.x+qq.y)/2, qq.y-qq.x, ylim=c(-1, 1),
     ylab="Private - Monopoly", xlab="mean")
abline(0,0)
title("Tukey's Mean Difference Plot")

#Kolmogorov-Smirnov test
ks.test(monopoly, private, alternative="two.sided")
ks.test(monopoly+0.5, private, alternative="two.sided")

#AUC
library(pROC)
ylab=factor(c(rep(0,length(monopoly)), rep(1,length(private))))
res = roc(ylab, c(monopoly, private))
res.ad = roc(ylab, c(monopoly+0.5, private))

par(mfrow=c(1,2))
plot(res, xlim=c(1,0), main="ROC between monopoly and private")
plot(res.ad, xlim=c(1,0), main="ROC between monopoly+0.5 and private")

#왼쪽그림의 AUC
auc(res) 
#오른쪽그림의 AUC
auc(res.ad)

#3.
x = c(rep(32,10),rep(33,23),rep(34,48),rep(35,80),rep(36,63),
      rep(37,65),rep(38,47),rep(39,33),rep(40,14),rep(42,6))

par(mfrow=c(2,2))
qqnorm(x, main='Normal prob plot')

n = length(x)
i = 1:n
q.exp.x = -log(1-(i-0.5)/n)          
plot(q.exp.x, x, main="Exponential prob plot")
plot(q.exp.x^(1/3), x^(1/3),main="Exponential; power=1/3")

q.weibull.x = log(q.exp.x)
plot(q.weibull.x, log(x), main="Weibull prob plot")

score.x = quantile(x, probs=c(0.25,0.75), names = FALSE,na.rm = TRUE)
score.y = qnorm(c(0.25,0.75))

slope = diff(score.x)/diff(score.y)
int = score.x[1] - slope * score.y[1]

slope
int

par(mfrow=c(1,1))
qqnorm(x, main='Normal prob plot')
abline(int, slope, lwd=2)

par(mfrow=c(1,1))
plot(q.weibull.x, log(x), main="Weibull prob plot")
abline(3.61148,0.04876, lwd=2)

(b=1/0.04876)
(a=exp(3.61148*(-b)))