1. 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?

# create data
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)
length(score)
## [1] 54
qqnorm(score, ylab="Score quantiles") ; qqline(score, col='red')
x=fivenum(score)
(pseudosigma = (x[4]-x[2])/1.34)
## [1] 2.238806
sd(score) # to compare
## [1] 2.092631
abline(x[3],pseudosigma,col="blue")
title(sub="intercept=13 (median); slope=2.239 (pseudo-sigma) blue line")

위 그림은 Score 자료에 대한 QQ-Plot이다. 빨간 선은 qqline()함수를 이용하여 그린 직선이고 파란 선은 수업 Script 파일을 참고하여 추정한 직선이다.

x <- fivenum(score)
x
## [1]  9 11 13 14 19

pseudo- sigma와 표준편차의 차이가 크지 않음을 먼저 확인하였고 이를 적용하여 추정하였다. 파란 선의 기울기는 pseudo- sigma를 대입하였고 절편은 median을 대입하였다. 두 선은 기울기는 거의 평행하지만 절편만 다름을 확인할 수 있다. 파란 선의 절편은 13, 빨간 선의 절편은 약 12.5로 추정되었다. 아마 절편의 값을 추정하는 부분에 있어 약간 다른 방식을 적용하는 듯하다.

하지만 결국 대부분의 점들이 두 직선위에 분포하고 있음을 확인할 수 있었고 score자료는 어느정도 정규분포를 따르고 있다고 판단하였다. 위 두 직선 중 파란선이 조금 더 잘 추정된 직선이라 판단하였고 이를 토대로 기울기와 절편을 판단한다면 2.239 / 13 이다.

hist(x = score, breaks = 10,freq=F)
lines(density(score),col="blue",lwd=3)
x=seq(-0,100,0.01)
y=dnorm(x,mean=13,sd=2.23) #normal density
lines(y~x,type='l',col="red",lwd=3)

y=dnorm(score,mean=13,sd=2.23)
stem(score,2)
## 
##   The decimal point is at the |
## 
##    9 | 0
##   10 | 0000
##   11 | 0000000000000
##   12 | 0000000
##   13 | 000000000
##   14 | 000000000
##   15 | 0000
##   16 | 00000
##   17 | 0
##   18 | 
##   19 | 0

Q-Q Plot을 이용하여 Score 데이터가 어느정도 정규분포를 따르고 있음을 확인하였고 실제 시각적으로 차이를 알아보고자 Score의 히스토그램과 줄기잎그림을 그렸다. 히스토그램에서 파란 곡선은 Score에 density를 함수를 적용하여 그린 선으로 이 데이터는 연속형이 아니기 때문에 시각적으로 참고하였다. 빨간 곡선은 실제로 평균:13, 표준편차:2.23인 정규분포를 생성하여 그렸다.

위 히스토그램과 줄기잎그림을 통해 수치적으롤 정확한 판단을 내리기는 힘들다. 데이터가 약간 왼쪽에 치우쳐 왜도가 보이긴 하지만 이러한 부분을 제외하고는 정규분포 곡선을 가지고 있다고 볼 수 있다. 실제 데이터가 정규분포를 따르고 있는지 수치상으로 확인할 수 있는 방법은 없을지 고민을 하였고 qnorm()함수를 이용하여 이를 확실하게 확인할 수 있었다.

quan=c(0,0.25,0.5,0.75,1)
qnorm(quan,13,2.23)
## [1]     -Inf 11.49589 13.00000 14.50411      Inf
fivenum(x)
## [1]   0  25  50  75 100

최소값과 최대값의 비교는 불가능하지만 제1분위수, 중간값, 제3분위수가 생성한 정규분포의 값과 비슷한 결과를 볼 수 있다. 또한, 이렇게 확인해보면서 제1분위수와 제3분위수 사이의 데이터들만 이용한다면 전체 데이터보다 더 확실하게 정규분포를 따를 것이라 예상해볼 수 있었다. 마지막으로 위의 정규분포의 평균과 표준편차는 위 직선의 추정식과 동일하게 13/2.239로 median 과 pseudo- sigma를 이용하였다.

다만 위 문제는 10개의 같은 값들로 이루어진 데이터라는 점과 데이터의 수가 많지 않다는 점이 판단하는 부분에 있어 어려움이 있었다.

2. Do the distributions of prices of two group follow the same distribution?

# create data
mono=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)
qqplot(mono,private,xlab="monopoly states",ylab="private-ownership states",
       main="Q-Q Plot");
abline(0,0.7)
abline(1,0.7)
abline(2,0.7)
abline(1.5,0.7)
abline(1.7,0.7)
abline(1.8,0.7)

#trial and error

Monopoly 자료와 private-ownership자료를 이용하여 QQ Plot을 얻었다. 코드를 통한 직선을 추정하기전에 직접 그림판을 이용하여 직선을 그어 기울기와 절편을 계산해보았다. 기울기는 Y축을 기준으로 0.2 증가할 때 X축이 얼마나 증가하는지 확인하여 구하였다. 위 그림의 경우 X값이 약4.1에서 4.38정도로 증가한다고 생각하였고 이를 통해 0.7 정도의 기울기를 얻었다. 절편의 경우 위 그림에서 구하는 것이 어려워 QQ Plot에서 기울기를 고정한 뒤 시행착오 과정을 통해 약 1.8의 값을 얻었다. abline()을 이용하여 Plot위에 추정한 직선(빨간색)과 코드를 이용한 직선(검정색) 2가지를 그렸고 그 결과는 아래와 같다.

abline(line(qqplot(mono,private,xlab="monopoly states"
                   ,ylab="private-ownership states",main="Q-Q Plot")));par(new=T)
abline(1.8,0.7,col="red")

직접 추정한 직선이 코드를 이용한 직선과 어느정도 맞아 떨어지는 것을 확인할 수 있다.

QQ Plot과 적절한 직선을 얻어 위처럼 시각적으로 표현하였다. QQ Plot을 확인해보면 X값이 4.2 부근에서는 직선위에 있다고 볼 수 있지만 4.2에서 멀어질수록 전체적인 데이터가 약간 직선에서 멀어지고 있다. 이는 분포의 모양이 완전히 같다고 할 수는 없지만 어느정도 비슷하다고 생각해볼 수 있다.

stem(mono)
## 
##   The decimal point is at the |
## 
##   3 | 8
##   4 | 00112222
##   4 | 566778
##   5 | 1
stem(private)
## 
##   The decimal point is 1 digit(s) to the left of the |
## 
##   42 | 99
##   44 | 045
##   46 | 55599
##   48 | 02555590555
##   50 | 0
##   52 | 0590
boxplot(mono,main="monopoly boxPlot")

boxplot(private,main="private boxPlot")

par(mfrow=c(2,1))
hist(mono,breaks = 15,probability = T,ylim =c(0,4),col="gray",xlim=c(3,6),
     xlab="monopoly states");par(new=T)
plot(density(mono),ylim =c(0,4),xlim=c(3,6),main="",col="red",lwd=2,xlab="")
hist(private,breaks = 15,probability = T,ylim =c(0,4),col="gray",xlim=c(3,6),
     xlab="private-ownership states");par(new=T)
plot(density(private),ylim =c(0,4),xlim=c(3,6),main="",col="blue",lwd=2,xlab="")

줄기잎그림, 히스토그램을 그려보았다. 줄기잎그림의 경우 비교하기가 힘들어 히스토그램만을 이용하여 비교를 진행하였다. 먼저 히스토그램을 그린 후 이를 바탕으로 density()를 통해 위와 같은 그림을 얻었다. 위 자료들은 연속형이 아니기 때문에 density를 나타낸 곡선은 시각적인 부분으로 참고하여야 할 것이다. 먼저 중간 값을 비교해보면 monopoly가 private-ownership보다 작을 것이다. 실제로 4.2, 4.85로 값이 차이가 난다.

또한, monopoly의 경우 자료가 정규분포처럼 보이지만 보다 퍼져 있는 느낌이 있다. 이는 데이터의 수가 적어서 나타날 수도 있겠다는 생각이 들었다. private-ownership은 가운데 부분은 정규분포 모양을 띄고 있지만 양 끝의 값들의 빈도도 꽤 높은 것을 확인할 수 있다. 마지막으로 X축의 길이를 동일하게 설정한 뒤 그렸기 때문에 데이터 분포 자체가 private-ownership이 좀 더 오른쪽으로 이동한 분포임을 쉽게 확인 가능하다.

마지막으로 상자그림을 그려보았는데 위에서 언급한 바와 같이 중간값 차이를 확인할 수 있었다. 또한, 히스토그램에서는 알아채지 못했던 특잇값이 private자료에 존재하고 있음을 확인할 수 있었다. 아마 이 특잇값으로 인해 아래의 QQ Plot에서 private의 그림이 양 끝으로 갈수록 직선과 멀어지는 결과가 나왔을 것이라 예상해볼 수 있다.

par(mfrow=c(1,2))
qqnorm(mono,main="monopoly normal QQ Plot");qqline((mono),col="red")
qqnorm(private,main="private normal QQ Plot");qqline((private),col="red")

위의 히스토그램을 보면서 두 분포 모두 정규분포의 모양이 보여 각 데이터의 QQ Plot을 그려보았다. 위에서 예상한 것처럼 monopoly 데이터의 점들이 어느정도 직선에 가까운 것을 확인할 수 있고 private자료는 가운데 부분은 확실하게 직선위에 분포하였지만 양 끝으로 갈수록 직선과 멀어지고 있음을 확인할 수 있었다.

data=c(mono,private)
qqnorm(data);qqline(data,col="red")

hist(data,freq=T,40,col="gray")

qqplot(mono,private,xlim=c(min(mono,private),max(mono,private)),
       ylim=c(min(mono,private),max(mono,private)))
title(main="QQ plot of mono, private")
abline(0,1,col="red")

첫번째 그림은 두 데이터를 합쳐 그린 히스토그램이다. 전 페이지에서 이야기한 것과 비슷하게 약간 거리가 있는 두 분포가 나타나고 있는 것을 확인할 수 있다. 두번째 그림은 qqnorm()을 이용하여 얻은 그림이다. 각자 다른 해석을 할 수 있겠지만 최대한 객관적으로 판단하고자 책 113p를 참고하여 여러 패턴과 비교해본 결과 혼합 정규분포로부터의 모의생성 자료에 대한 정규 확률 플롯과 가장 유사하다고 생각하였다. 물론 두 데이터가 정규분포라는 강력한 증거는 없지만 위의 QQ Plot에서 어느정도 정규분포 가정이 가능하다 생각하여 두 데이터를 합쳐 위와 같은 그림을 그렸다. 위 Plot을 통해 두 분포의 모수가 차이가 있다는 판단을 내리게 되었다.

(qqplot(mono,private))
## $x
##  [1] 3.80 4.00 4.00 4.10 4.11 4.15 4.19 4.20 4.20 4.50 4.55 4.55 4.65 4.74 4.75
## [16] 5.05
## 
## $y
##  [1] 4.290000 4.430000 4.543333 4.750000 4.750000 4.790000 4.800000 4.840000
##  [9] 4.850000 4.850000 4.896667 4.950000 4.950000 5.166667 5.263333 5.300000
qq.x <- qqplot(mono,private)$x

qq.y <- qqplot(mono,private)$y

x11()
plot((qq.x+qq.y)/2, qq.y-qq.x, main="Tukey mean difference plot",ylim=c(-1,1),
     ylab="private-monopoly", xlab="mean")
abline(0,0)

마지막으로 Tukey’s Mean-Difference plot을 이용하여 차이를 확인하였다.

모든 점들의 값이 0 이상에 분포하고 있다. 즉, 차이가 있다는 것이다. 위의 과정을 토대로 결론을 내리자면 두 분포가 정규분포와 비슷한 모양을 띄고 있으며 디테일한 부분에서 약간씩 차이가 있다는 것을 확인하였다. 또한 두 분포의 중간값이 어느정도 차이가 났기 때문에 분포가 나타나는 위치가 달랐다. 마지막으로 monopoly states의 데이터가 private-ownership states에 비해 적기 때문에 두 분포의 차이를 확인하는데 어려움이 있었다. 히스토그램을 살펴보았을 때도 하나의 값이 전체 데이터에 큰 영향을 미치고 있어 private-ownership states에 비해 뚜렷한 분포가 나타나지 않아 판단에 모호함이 존재하였다.

3. 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?

# create data
x=c(32,33,34,35,36,37,38,39,40,42)
frequency=c(10,33,81,161,224,289,336,369,383,389)
x1=c()
for(i in 1:length(x)){
        x1[i]=frequency[i]-frequency[i-1]        
}
x1[1]=10
strength=c()
for(i in 1:length(x)){
  strength=c(strength,(rep(x[i],x1[i]))) 
}

주어진 데이터를 각 빈도수로 바꾸어 얻은 자료의 히스토그램, boxplot은 아래와 같다.

hist(strength,breaks = 10,probability = T,xlim=c(25,50),col="gray")
lines(density(strength),col="blue",lwd=2)

boxplot(strength,main="yield strength of a Bofores steel")

어느정도 정규분포 형태를 보이는 듯하며 여러 QQ Plot을 통해 어떤 분포를 따르는지 확인해 보자.

n <- length(strength)
i <- 1:n
q.exp.strength <- -log(1-(i-0.5)/n)             # exponential quantile
plot(q.exp.strength, strength, main="Exponential prob plot")

line(qqplot(q.exp.strength, strength))
## 
## Call:
## line(qqplot(q.exp.strength, strength))
## 
## Coefficients:
## [1]  33.92   2.49
abline(33.92,2.49) #qqline
abline(32,3,col="red")
abline(35,2,col="red")
abline(34.5,2,col="red")

자료의 지수분포 확률 Plot은 위와 같다. 검은 선은 qqline()을 이용하여 그은 선이고 빨간 선은 여러 시도 끝에 얻은 가장 적절하다고 생각되는 선이다. 하지만 그림의 결과와 같이 많은 선을 그어보아도 적합한 선을 찾을 수 없었다. Plot 자체도 시각적으로 곡선을 형태를 띄고 있는 듯하여 지수분포는 위의 자료는 지수분포를 따르지 않는다는 판단을 하였다.

plot(q.exp.strength^(1/3), strength^(1/3),main="Exponential prob plot",
     sub="re-expression; power=1/3")

line(qqplot(q.exp.strength^(1/3), strength^(1/3)))
## 
## Call:
## line(qqplot(q.exp.strength^(1/3), strength^(1/3)))
## 
## Coefficients:
## [1]  3.1369  0.1893
abline(3.14,0.19)
abline(3.135,0.19,col="red")

1/3승을 적용하여 재표현한 지수확률 분포 Plot은 위와 같다. 검은 선은 qqline()을 이용하여 그은 선이고 빨간 선은 여러 시도 끝에 얻은 가장 적절하다고 생각되는 선이다. 대부분의 점들이 직선위에 분포하고 있다. 재표현하기 전의 Plot에 비하면 훨씬 더 적합하다고 판단되며 위 데이터가 어떤 분포인지 판단할 때 충분히 고려해 볼만하다고 생각하였다.

q.weibull.strength <- log(q.exp.strength)
x11()
plot(q.weibull.strength, log(strength), main="Weibull prob plot")

line(qqplot(q.weibull.strength,log(strength)))
## 
## Call:
## line(qqplot(q.weibull.strength, log(strength)))
## 
## Coefficients:
## [1]  3.61148  0.04876
abline(3.611,0.0487)
abline(3.605,0.04,col="red")

와이블 분포 확률 Plot은 위와 같다. 검은 선은 qqline()을 이용하여 그은 선이고 빨간 선은 여러 시도 끝에 얻은 가장 적절하다고 생각되는 선이다. 지수분포 확률 Plot을 뒤집어 놓은 형태처럼 보이는 듯하지만 간격이 조금 다르게 나타나고 있다. 이 Plot 또한 어떠한 직선에도 모든 점들이 분포하지 않고 있으며 점들이 지수 분포 확률 Plot과 비슷하게 약간 곡선의 형태를 가지고 있는 것처럼 보인다. 이러한 이유로 이 데이터의 분포는 와이블 분포를 따를 가능성이 낮다고 판단하였다. 현재는 점들이 직선에 분포하고 있지 않지만 만약 분포하는 데이터를 얻게 된다면 형태(shape)를 나타내는 a와 척도(scale)을 나타내는 b의 값은 각각 검은 선을 기준으로 a=e^(intercept*(-b)),b=1/slope 이다.

이를 계산하면 다음과 같은 값을 얻을 수 있다.

(b=1/0.0487)
## [1] 20.53388
(a=exp(3.611*(-b)))
## [1] 6.280591e-33
qqnorm(strength) #normal check
qqline(strength)
x <- fivenum(strength)
x
## [1] 32 35 36 38 42
(pseudosigma = (x[4]-x[2])/1.34)
## [1] 2.238806
median(strength)
## [1] 36
abline(36,2.1,col="red")

마지막으로 자료에 대한 정규확률 Plot이다. 검은 선은 qqline()을 이용하여 그은 선이고 빨간 선은 여러 시도 끝에 얻은 가장 적절하다고 생각되는 선이다. 이전의 와이블 분포와 지수 분포 확률 Plot에 비해 대부분의 점들이 선에 위치하고 있다. 여기서는 빨간 선이 더 적합하다고 판단하여 이를 기준선으로 삼았다. 빨간 선의 기울기는 2.1, 절편은 36이다. 기울기는 pseudo-sigma를 이용하여 어느정도 기준을 삼았으며 절편은 중간값을 사용하였다.

지금까지 데이터가 어느 분포를 가장 잘 따르고 있는지 확인해보았다. 예상되는 후보로는 1/3승을 적용하여 재표현한 지수확률 분포와 정규확률 분포이다. 이중 반드시 하나의 분포를 택해야 한다면 정규확률 분포를 택할 것이다. 그 이유는 처음 히스토그램을 통해 데이터가 어느 정도 정규분포를 따르고 있음을 확인하였고 기울기와 절편을 추정하는 것이 더 간편하다고 생각되기 때문이다.

그뿐 아니라 정규분포가 익숙한 개념이기 때문에 쉽게 분석할 수 있을 것이다. 또한 1번에서처럼 qnorm()을 이용하여 제1사분위수, 중간값, 제3사분위수가 실제 데이터와 어느정도 맞아 떨어졌기 때문이다.