Re-expression
Exploratory Data Analysis
Choi Jong Moon
21 April 2020
1. DISTRESS.DAT는 Severe Idiopathic Respiratory Distress Syndrome이 있는 신생아의 태어날 때 몸무게(Kg)이다. (*표는 죽은 아이이다. 자료파일은 ASCII 형식이다.)
(자료가 5개 칸으로 되어 있는 것은 특별한 의미 없다. 5개 집단 아님)
(사망, 생존 여부를 변수로 만들어 편집 한 후 R로 읽는 것이 편하다.)
가. 생존 여부에 관계없이 전체 몸무게의 자료로 상자그림과 줄기그림을 그리고 대칭화 재표현을 한 후의 그림들과 비교하여라.
“DISTRESS.dat” 파일을 불러와 data라는 변수로 저장한다. 자료는 5개 칸으로 되어있지만 이것이 집단을 의미하는 것은 아니므로 이를 벡터화한다. 생존 여부에 관계없이 데이터를 분석하기 위해 * 표시를 없애준다. 자료를 분석하기 위해 모든 값을 numeric으로 바꿔준다. 전처리한 최종 데이터는 data2라는 변수에 저장했으며 그 결과는 다음과 같다.
data = read.table("DISTRESS.dat", header=FALSE)
data = as.vector(as.matrix(data))
library(stringr)
data2 = str_replace_all(data, "\\*", "")
data2 = as.numeric(data2)
data2## [1] 1.050 1.175 1.230 1.310 1.500 1.600 1.720 1.750 1.770 2.275 2.500 1.030
## [13] 1.100 1.185 1.225 1.262 1.295 1.300 1.550 1.820 1.890 1.940 2.200 2.270
## [25] 2.440 2.560 2.730 1.130 1.575 1.680 1.760 1.930 2.015 2.090 2.600 2.700
## [37] 2.950 3.160 3.400 3.640 2.830 1.410 1.715 1.720 2.040 2.200 2.400 2.550
## [49] 2.570 3.005대칭화 재표현이란 f(x)=x^p에 대해 자료를 대칭적으로 만드는 최적의 p를 찾는 작업이다. Taylor Expansion을 2nd order term까지 정리해 조건식 Hu-M=M-Hl으로 풀어주면 다음과 같은 p를 도출할 수 있다.
p = 1 - (2 * M * ( Hu-M + Hl-M ) / ( (Hl-M)^2 + (Hu-M)^2 ))
M은 Median, Hu는 Upper Hinge, Hl은 Lower Hinge를 의미하며 조건식에서의 각각의 값은 x^p를 적용한 자료에 대한 정보이고, p를 계산할 때의 각각의 값은 원자료에 대한 정보이다.
p를 구하기 위해 다섯수치요약을 수행한다. Upper Hinge는 2.5, Lower Hinge는 1.41이다.
fivenum(data2)## [1] 1.030 1.410 1.855 2.500 3.640p계산을 수행해주는 re_exp() 함수를 정의한다. 자료를 벡터형태로 입력하면 최적의 p를 도출해준다. 함수에 신생아 몸무게의 자료를 넣었더니 -0.2083707이라는 결과가 나왔다. 계산된 p를 그대로 가져다쓰는 것보다 이를 표준화된 숫자로 바꿔주는 것이 분석에 용이하다. 따라서 p=0이라고 가정하여 신생아 몸무게에 대한 재표현으로 로그변환을 채택한다
re_exp = function(x){
Hl <- fivenum(x)[2]
M <- fivenum(x)[3]
Hu <- fivenum(x)[4]
return(1-2*M*(Hu-M+Hl-M)/((Hl-M)^2+(Hu-M)^2))
}
re_exp(data2)## [1] -0.2083707이제 원자료와 로그변환한 자료에 대해 각각 줄기그림과 상자그림을 그려보자. 먼저 각각에 대한 줄기그림은 다음과 같다.
# 신생아 몸무게 줄기잎그림
stem(data2)##
## The decimal point is at the |
##
## 1 | 0111222233334
## 1 | 566677778888999
## 2 | 001223344
## 2 | 56666778
## 3 | 0024
## 3 | 6# 로그 신생아 몸무게 줄기잎그림
stem(log(data2))##
## The decimal point is 1 digit(s) to the left of the |
##
## 0 | 350267
## 2 | 0136674
## 4 | 14572444677
## 6 | 046601499
## 8 | 2289244469
## 10 | 04805
## 12 | 29두 개의 줄기잎그림을 살펴보면 먼저 가장 큰 차이는 줄기에 있다. 자료에 로그를 취해주었기 때문에 자료의 기준이 되는 줄기의 단위가 매우 달라졌다. 따라서 값에 대한 비교는 불가능하며 자료의 전체적인 분포 모양만을 비교하는 것이 적절하다. 각각은 모두 하나의 큰 군집으로 이루어져있다. 별다른 특이값도 보이지 않는다. 역시 가장 큰 차이는 대칭성 여부이다. 원자료는 오른쪽으로 길게 늘어진 형태의 분포를 띤다. 하지만 로그를 취해준 자료는 종모양의 대칭적 분포를 띰을 확인할 수 있다. 대칭화 재표현이 적절하게 이루어졌다. 대칭화를 통해 mean과 median이 비슷해졌음을 유추할 수 있다. 또한 분포가 종모양을 따르므로 근사적으로 정규분포에 가까운 모양이다.
## [1] "원자료 평균: 1.97494"## [1] "원자료 중앙값: 1.855"## [1] "로그 변환 자료 평균: 0.625626484857429"## [1] "로그 변환 자료 중앙값: 0.617706665080128"
이번에는 두 자료에 대해 상자그림을 그려보자. 결과는 다음과 같다. 
두 그림의 가장 큰 차이는 수염의 길이이다. 원자료의 상자그림은 위아래로 뻗은 수염의 길이에 차이가 크다. 하지만 로그변환 후의 상자그림은 위아래의 수염의 길이가 한눈에도 비슷해보인다. 즉 대칭화가 올바르게 이루어졌다는 것이다. 또한 중앙값을 나타내는 검은색 선이 상자의 한가운데로 이동했다. 검은색 선이 상자의 중앙에서 멀어질수록 왜도가 심해진다. 로그변환 후의 상자그림에서는 중앙값이 상자의 거의 한가운데 위치하므로 역시 대칭화가 잘 이루어졌음을 확인할 수 있다. 실제로 아래의 식을 통해 왜도를 확인해보면 로그변환 후의 데이터의 왜도는 거의 0에 가까움을 확인할 수 있다.
(EDA적으로 왜도를 계산하는 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(data2)## [1] 0.1834862# 로그변환 자료 왜도
skewness(log(data2))## [1] 0.0427223나. 사망과 생존 집단으로 나누어 상자 그림을 그리고 산포(spread)가 같도록 하는 재표현을 찾아 상자 그림을 그리고 두 집단을 비교하여라. 위의 병이 있는 신생아 중 몸무게가 몇 Kg이하면 사망할 가능성이 많다고 예측할 수 있겠는가?
이제는 사망과 생존 집단으로 나누어 자료를 비교해보자. 먼저 이를 위해 새로 데이터프레임을 구성했다. 기존 data자료에서 27번째까지의 값은 사망, 28번째부터의 값은 생존이다. 각각의 label을 생성하고 data2와 합쳐 df라는 새로운 데이터프레임을 생성했다. df는 아래와 같다.
label = c(rep("사망", 27), rep("생존", length(data)-27))
df = data.frame(data2, label)
colnames(df) = c("몸무게", "생존여부" )
head(df)각각의 자료에 대해 대칭화 재표현을 해주는 최적의 p를 찾기 위해re_exp()를 적용한다.
# 생존 신생아의 최적의 p
re_exp(df[df$생존여부=="생존",]$몸무게)## [1] 0.1296567# 사망 신생아의 최적의 p
re_exp(df[df$생존여부=="사망",]$몸무게)## [1] -0.07216304각각은 0.13과 -0.07 정도가 나왔으므로 모두 표준화된 숫자인 0으로 가정하여 로그를 취하는 것이 바람직하다. 하지만 이번에는 산포를 같게 하는 재표현을 찾아보도록 하자.
산포를 같게 하는 재표현은 log(H-spread) = k + b * log(M)을 그래프로 그려 기울기 b를 찾는 일이다. 회귀식을 통해 기울기를 계산해도 되고 눈대중으로 기울기를 확인해도 된다. 마지막으로는 1-b를 계산하여 최적의 p를 찾는다. 다음은 회귀식의 결과와 그래프이다. 이를 통해 b를 계산할 수 있다.
##
## Call:
## lm(formula = log(Hspr) ~ log(M))
##
## Coefficients:
## (Intercept) log(M)
## -0.5157 0.6854
최종적으로 확인한 b=0.6854이므로 p=1-b=0.3146이다. 이를 반올림하면 0.5이다. 0.5는 Square Root 변환이다. 실제로 이 방법을 통해 산포(H-spread)의 차이가 줄었는지를 확인해보자.
HsprA <- fivenum(df[df$생존여부=="사망",]$몸무게)[4] - fivenum(df[df$생존여부=="사망",]$몸무게)[2]
HsprB <- fivenum(df[df$생존여부=="생존",]$몸무게)[4] - fivenum(df[df$생존여부=="생존",]$몸무게)[2]
HsprA - HsprB## [1] -0.201HsprC <- fivenum(df[df$생존여부=="사망",]$몸무게^0.5)[4] - fivenum(df[df$생존여부=="사망",]$몸무게^0.5)[2]
HsprD <- fivenum(df[df$생존여부=="생존",]$몸무게^0.5)[4] - fivenum(df[df$생존여부=="생존",]$몸무게^0.5)[2]
HsprC - HsprD## [1] -0.02182665원자료의 생존여부에 따른 H-spread의 차이는 -0.201이다. 루트 변환 후의 생존여부에 따른 H-spread의 차이는 -0.02182665이다. 산포의 차이가 많이 줄어들었음을 확인할 수 있다.
각각의 자료에 Square Root를 취해주어 상자그림을 그려보면 다음과 같다. 
정확한 비교를 위해 notch를 포함한 상자그림을 그렸다. 상자그림에서 notch는 각각의 상자에서 안쪽으로 꺾이는 범위를 통해 확인할 수 있다. notch는 EDA에서 사용하는 robust한 신뢰구간으로 중앙값에 대한 95% 신뢰구간을 계산한다. 비교하는 대상의 notch의 범위가 중첩하지 않는다면 대략 95%의 유의수준에서 두 자료는 유의미한 차이가 있다고 결론내릴 수 있다.
자료를 비교하기 위해서는 산포를 같게 해야한다. 위에서는 산포를 같게 하는 방법으로 p=0.5를 찾았고 오른쪽 상자그림이 Square Root를 취해준 상자그림이다. 대칭화의 결과로 중앙값이 상자의 중앙으로 이동했으며 상자의 수염의 길이가 비슷해졌다. 이를 기준으로 결과를 분석해보자. 결과를 보면 사망한 신생아의 로그몸무게 중앙값은 1.25정도이다. 생존한 신생아의 루트몸무게 중앙값은 1.5정도이다. 약 95%의 유의 수준에서 사망한 신생아의 몸무게 중앙값은 생존한 신생아의 몸무게 중앙값보다 낮다고 말할 수 있다. 몸무게가 낮을수록 사망할 가능성이 높다는 것이다. 그렇다면 그 기준은 몇kg이라고 할 수 있을까? 이 경우에는 각각의 notch가 1.35정도를 기준으로 위아래로 나뉘기 때문에 그 기준을 1.35로 정함이 바람직하다. 루트몸무게가 1.35이므로 원자료로 환산하면 1.35^2=1.8225이다. 즉, 태어날 때의 몸무게가 1.8225kg 미만이라면 사망할 가능성이 높다고 말할 수 있다.
2. R의 airquality 자료에서 Ozone 자료를 대칭화하기 위한 변환공식을 이용하여 변환 전과 변화 후를 비교하여라.
airquality는 다음의 정보를 갖는 데이터프레임이다. 데이터셋에 대한 설명은 ?airquality라는 명령어를 통해 확인할 수 있다.
New York Air Quality Measurements
- Description : Daily air quality measurements in New York, May to September 1973.
- Format : A data frame with 153 observations on 6 variables.
- Columns
[,1] Ozone : numeric Ozone (ppb)
[,2] Solar.R : numeric Solar R (lang)
[,3] Wind : numeric Wind (mph)
[,4] Temp : numeric Temperature (degrees F)
[,5] Month : numeric Month (1–12)
[,6] Day : numeric Day of month (1–31)
대칭화 재표현이란 f(x)=x^p에 대해 자료를 대칭적으로 만드는 최적의 p를 찾는 작업이다. Taylor Expansion을 2nd order term까지 정리해 조건식 Hu-M=M-Hl으로 풀어주면 다음과 같은 p를 도출할 수 있다.
p = 1 - (2 * M * ( Hu-M + Hl-M ) / ( (Hl-M)^2 + (Hu-M)^2 ))
M은 Median, Hu는 Upper Hinge, Hl은 Lower Hinge를 의미하며 조건식에서의 각각의 값은 x^p를 적용한 자료에 대한 정보이고, p를 계산할 때의 각각의 값은 원자료에 대한 정보이다.
p를 구하기 위해 다섯수치요약을 수행한다. Upper Hinge는 63.5, Lower Hinge는 18이다.
fivenum(airquality$Ozone, na.rm=T)## [1] 1.0 18.0 31.5 63.5 168.0위에서 정의한 p계산을 수행해주는 re_exp() 함수를 사용한다. 자료를 벡터형태로 입력하면 최적의 p를 도출해준다. 함수에 Ozone의 자료를 넣었더니 0.03378238이라는 결과가 나왔다. 계산된 p를 그대로 가져다쓰는 것보다 이를 표준화된 숫자로 바꿔주는 것이 분석에 용이하다. 따라서 p=0이라고 가정하여 신생아 몸무게에 대한 재표현으로 로그변환을 채택한다
re_exp(airquality$Ozone)## [1] 0.03378238원자료와 로그 변환 자료에 대해 각각 상자그림을 그린다. 
수치적 비교를 위해 각각의 평균, 중앙값과 왜도를 구한다. 왜도는 위에서 정의한 EDA적 skewness() 함수를 사용한다.
## [1] "원자료 평균: 42.1293103448276"## [1] "원자료 중앙값: 31.5"## [1] "로그 변환 자료 평균: 3.41851510081201"## [1] "로그 변환 자료 중앙값: 3.44986155364244"# 원자료 왜도
skewness(airquality$Ozone)## [1] 0.4065934# 로그 변환 자료 왜도
skewness(log(airquality$Ozone))## [1] 0.1123698위의 내용을 종합하면 대칭화가 잘 이루어졌음을 확인할 수 있다. 변환 후, 중앙값이 상자의 중앙으로 이동했으며 상자의 수염의 길이가 비슷하게 조정되었다. 이를 수치로 확인할 수 있는데 원자료의 평균과 중앙값의 차이보다 로그 변환 자료의 평균과 중앙값의 차이가 현격히 작아졌음을 확인할 수 있다.(물론 값의 크기에 대해 스케일링을 해주어도 후자가 더 작다.) 또한 Hinge를 통해 계산한 EDA적 왜도값 역시 원자료는 0.41, 로그 변환 후 자료는 0.11로 후자가 훨씬 더 0에 가깝다. 대칭화가 올바르게 이루어졌다는 것을 의미한다. 다만 상자그림에서 주목할 점은 높은 값에 존재했던 두 개의 특이값이 사라지고 낮은 값 쪽에서 하나의 특이값이 새로 발생했다는 점이다. 로그는 큰 값의 범위를 작게 만들어주기 때문에 기존의 특이값을 분석하여 이를 변환 후에도 특이값으로 처리해야 하는지를 확인해야 한다.
3. 2016년의 삼성전자 주식 가격을 인터넷에서 찾아 줄기잎 전시와 상자그림을 그려 대칭성을 점검하고 필요하다면 대칭화하는 변환을 시행착오방법으로 찾아라. 어떤 시행 착오로 찾아 갔는지 중간 과정을 전부 보여라. 변환 후의 줄기잎 전시와 상자그림을 그려 전, 후를 비교하여라. 변환 전, 후의 skewness 값을 비교하여라.
아래는 네이버 금융 웹페이지에서 삼성전자의 주가 정보를 크롤링하는 코드이다. 크롤링을 통해 2016년의 삼성전자 주가 데이터를 확보할 수 있었다.
#Load Data
library(rvest)
library(httr)
library(stringr)
ticker=list()
for (i in 1:150) {
url=paste0("https://finance.naver.com/item/sise_day.nhn?code=005930&page=",i)
down = GET(url)
Sys.setlocale("LC_ALL", "English")
table = read_html(down, encoding = "EUC-KR") %>%
html_table(., fill = TRUE)
table = table[[1]]
table = table[-c(1,7,8,9,15),]
Sys.setlocale("LC_ALL", "Korean")
ticker[[i]] = table
}
ticker = do.call(rbind, ticker)
rownames(ticker) = NULL
ticker_2016 = ticker[grep("^2016", ticker$날짜),]
rownames(ticker_2016) = NULL최종적으로 얻은 데이터는 다음과 같다.
head(ticker_2016)우리가 원하는 정보는 종가이므로 이를 numeric으로 전처리하여 price라는 변수에 저장한다.
par(mfrow=c(1,1))
price = as.numeric(str_replace_all(ticker_2016$종가, ",",""))
print(head(price))## [1] 1802000 1788000 1799000 1798000 1782000 1809000
재표현을 수행하기 전에 줄기잎그림과 상자그림을 통해 자료를 살펴보자.
stem(price)##
## The decimal point is 5 digit(s) to the right of the |
##
## 11 | 33333444
## 11 | 555555566666677777788888899999
## 12 | 0112223
## 12 | 555555566666677777777778888888889999999999
## 13 | 00000000001113
## 13 | 7788
## 14 | 000001111233333
## 14 | 55566777789
## 15 | 001222233333444444
## 15 | 5555666667777777777778899999999
## 16 | 00000011111222222334444444444
## 16 | 55555557888899
## 17 | 123
## 17 | 5555677888999
## 18 | 0000111
줄기 그림을 통해 크게 3개의 군집이 존재함을 확인할 수 있었다. 상자 그림에서는 중앙값이 제1사분위수보다는 제3사분위수에 조금 더 가깝다. 자료의 분포는 꽤 대칭적이지만 왼쪽으로 살짝 늘어진 형태를 보이고 있다.
그동안의 문제에서는 re_exp()함수를 활용한 대칭화 재표현, 산포를 줄이는 재표현의 방법으로 재표현 방법을 채택했다. 이번 문제에서는 시행착오적 방법으로 적절한 재표현 방법을 찾아보자. 시행착오적 방법이란 f(x) = x^p에 다양한 p를 직접 넣어가며 가장 적절한 p를 찾는 방법이다.
자료를 대칭화하고 싶으므로 왜도를 기준으로 p를 찾아보자. 다양한 p에 대해 위에서 정의한 skewness() 함수를 적용하여 왜도값을 산출한다. -5부터 5까지 0.1의 간격으로 p의 후보들을 생성한다. 그 후에 자료에 p승을 하여 각각의 모든 왜도값을 산출한다. 또한 우리는 왜도의 방향이 아닌 크기를 확인하고 싶으므로 이에 대해 모두 절대값을 취해준다.
p = seq(-5, 5, by=0.1)
skewness_p = c()
for (i in p){
n = skewness(price^i)
skewness_p = c(skewness_p, n)
}
# p=0에 대해 log값을 넣어준다.
skewness_p[51] = skewness(log(price))
# 절대값을 취한다.
skewness_p = abs(skewness_p)이를 다음과 같이 그래프로 표현할 수 있다. 
p[which(skewness_p == min(skewness_p))]## [1] 3.4결과적으로 가장 작은 왜도를 산출하는 p는 3.4이다. 이를 적용하여 재표현된 줄기잎그림과 상자그림을 그려보자.
#원자료
stem(price)##
## The decimal point is 5 digit(s) to the right of the |
##
## 11 | 33333444
## 11 | 555555566666677777788888899999
## 12 | 0112223
## 12 | 555555566666677777777778888888889999999999
## 13 | 00000000001113
## 13 | 7788
## 14 | 000001111233333
## 14 | 55566777789
## 15 | 001222233333444444
## 15 | 5555666667777777777778899999999
## 16 | 00000011111222222334444444444
## 16 | 55555557888899
## 17 | 123
## 17 | 5555677888999
## 18 | 0000111#재표현자료
stem(price^3.4)##
## The decimal point is 20 digit(s) to the right of the |
##
## 3 | 88888999
## 4 | 00000111122223333333334445556667899
## 5 | 00333344445556666666667777888888889999
## 6 | 00000001111111112337
## 7 | 235588899
## 8 | 0011344555889
## 9 | 12222357
## 10 | 00134446677789999
## 11 | 000123334456666666667788
## 12 | 0111122233344466788899
## 13 | 0001155555556666888889
## 14 | 256679
## 15 | 048
## 16 | 17889
## 17 | 146789
## 18 | 123455678
## 19 | 03.4승으로 재표현한 줄기잎그림의 가장 큰 특징은 오른쪽으로 긴 꼬리가 생겼다는 점이다. 왜도에서는 거의 0에 수렴하는 자료가 실제로는 이렇게 꼬리가 긴 이유는, EDA의 왜도는 H-spread내에서의 자료만으로 계산하기 때문이다. 따라서 그 바깥의 자료에 대해서는 대칭이라고 할 수 없다. 이는 상자그림을 통해서도 확인 가능하다. 
상자그림을 한눈에 보면 오히려 원자료가 더 대칭적인 것처럼 보인다. 재표현된 자료는 높은 쪽의 꼬리가 너무 길어졌다. 자료에 3.4승을 해주었기 때문에 큰 값이 더 커진 것이다. 다만, 상자 내부를 살펴보면 중앙값이 거의 완벽하게 상자 중앙에 위치함을 확인할 수 있다. 상자 내부에서는 거의 완벽한 대칭이 이루어졌다.
왜도를 기준으로 하는 p의 선정에는 위와 같은 단점이 발생하므로 다른 방법을 찾아보자. 눈으로 확인하는 것이다. 대표적인 변환인 p=0(로그), p=0.5, p=-1, p=-0.5, p=2에 대해 상자그림을 그려보자. 그림에 각 분포의 왜도값도 표시했다. 
로그와 루트 변환은 수염의 길이는 조정이 되었지만 중앙값이 상자의 중앙에서 조금 더 멀어지는 단점이 발생한다. 역수와 루트역수의 변환도 적절한 대칭이라고는 보기 힘들다. 오히려 원자료의 왜도가 0.13으로 이 중 가장 낮다. p=2에서는 중앙값은 상자의 중앙으로 이동했지만 위쪽 수염이 길어진다는 단점이 있다. 즉 상자내부의 최빈값들의 정보가 중요하다면 3.4승의 재표현을 통해 대칭적인 분포를 만들어야 한다. 하지만 이미 원자료의 중앙값도 꽤 상자의 중앙에 위치하기 때문에 굳이 수염의 비대칭을 감수할 필요는 없다. 따라서 이 자료는 원자료 그대로 두는 것이 가장 적절하다.
R Script
보고서에 쓰인 모든 코드를 아래에 정리해두었다.
#1.
data = read.table("DISTRESS.dat", header=FALSE)
data = as.vector(as.matrix(data))
#1.가.
library(stringr)
data2 = str_replace_all(data, "\\*", "")
data2 = as.numeric(data2)
data2
fivenum(data2)
re_exp = function(x){
Hl <- fivenum(x)[2]
M <- fivenum(x)[3]
Hu <- fivenum(x)[4]
return(1-2*M*(Hu-M+Hl-M)/((Hl-M)^2+(Hu-M)^2))
}
re_exp(data2)
stem(data2)
stem(log(data2))
paste("원자료 평균:", as.character(mean(data2)))
paste("원자료 중앙값:", as.character(median(data2)))
paste("로그 변환 자료 평균:", as.character(mean(log(data2))))
paste("로그 변환 자료 중앙값:", as.character(median(log(data2))))
par(mfrow=c(1,2))
boxplot(data2, col="papayawhip", ylab="몸무게")
title("Boxplot of Weights")
boxplot(log(data2), col="papayawhip", ylab="로그몸무게")
title("Boxplot of Log Weights")
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(data2)
# 로그변환 자료 왜도
skewness(log(data2))
#1.나
label = c(rep("사망", 27), rep("생존", length(data)-27))
df = data.frame(data2, label)
colnames(df) = c("몸무게", "생존여부" )
head(df)
re_exp(df[df$생존여부=="생존",]$몸무게)
re_exp(df[df$생존여부=="사망",]$몸무게)
HsprA <- fivenum(df[df$생존여부=="사망",]$몸무게)[4] - fivenum(df[df$생존여부=="사망",]$몸무게)[2]
HsprB <- fivenum(df[df$생존여부=="생존",]$몸무게)[4] - fivenum(df[df$생존여부=="생존",]$몸무게)[2]
M <- c(median(df[df$생존여부=="사망",]$몸무게), median(df[df$생존여부=="생존",]$몸무게))
Hspr = c(HsprA, HsprB)
(RegrLine <- lm(log(Hspr) ~ log(M))) # ( ) is needed to print coeff's
plot(log(M), log(Hspr), main="Spread vs. Level plot")
abline(coef(RegrLine))
HsprA <- fivenum(df[df$생존여부=="사망",]$몸무게)[4] - fivenum(df[df$생존여부=="사망",]$몸무게)[2]
HsprB <- fivenum(df[df$생존여부=="생존",]$몸무게)[4] - fivenum(df[df$생존여부=="생존",]$몸무게)[2]
HsprA - HsprB
HsprC <- fivenum(df[df$생존여부=="사망",]$몸무게^0.5)[4] - fivenum(df[df$생존여부=="사망",]$몸무게^0.5)[2]
HsprD <- fivenum(df[df$생존여부=="생존",]$몸무게^0.5)[4] - fivenum(df[df$생존여부=="생존",]$몸무게^0.5)[2]
HsprC - HsprD
par(mfrow=c(1,2))
boxplot(몸무게 ~ 생존여부, data=df, col="papayawhip", notch=T)
title("Boxplot of Weights")
df$루트몸무게 = sqrt(df$몸무게)
boxplot(루트몸무게 ~ 생존여부, data=df, col="papayawhip", notch=T)
title("Boxplot of Square root Weights")
#2.
fivenum(airquality$Ozone, na.rm=T)
re_exp(airquality$Ozone)
par(mfrow=c(1,2))
boxplot(airquality$Ozone, na.rm=T, col="papayawhip")
title("Boxplot of Ozone")
boxplot(log(airquality$Ozone), na.rm=T, col="papayawhip")
title("Boxplot of Log Ozone")
paste("원자료 평균:", as.character(mean(airquality$Ozone, na.rm=TRUE)))
paste("원자료 중앙값:", as.character(median(airquality$Ozone, na.rm=TRUE)))
paste("로그 변환 자료 평균:", as.character(mean(log(airquality$Ozone), na.rm=TRUE)))
paste("로그 변환 자료 중앙값:", as.character(median(log(airquality$Ozone), na.rm=TRUE)))
# 원자료 왜도
skewness(airquality$Ozone)
# 로그 변환 자료 왜도
skewness(log(airquality$Ozone))
#3.
#Load Data
library(rvest)
library(httr)
library(stringr)
ticker=list()
for (i in 1:150) {
url=paste0("https://finance.naver.com/item/sise_day.nhn?code=005930&page=",i)
down = GET(url)
Sys.setlocale("LC_ALL", "English")
table = read_html(down, encoding = "EUC-KR") %>%
html_table(., fill = TRUE)
table = table[[1]]
table = table[-c(1,7,8,9,15),]
Sys.setlocale("LC_ALL", "Korean")
ticker[[i]] = table
}
ticker = do.call(rbind, ticker)
rownames(ticker) = NULL
ticker_2016 = ticker[grep("^2016", ticker$날짜),]
rownames(ticker_2016) = NULL
par(mfrow=c(1,1))
price = as.numeric(str_replace_all(ticker_2016$종가, ",",""))
print(head(price))
#stem and boxplot
par(mfrow=c(1,1))
stem(price)
boxplot(price, col="papayawhip", ylim=c(1000000, 2000000), ylab="price(won)")
title("Boxplot of Samsung Electronics' stock price(2016)")
p = seq(-5, 5, by=0.1)
skewness_p = c()
for (i in p){
n = skewness(price^i)
skewness_p = c(skewness_p, n)
}
# p=0에 대해 log값을 넣어준다.
skewness_p[51] = skewness(log(price))
# 절대값을 취한다.
skewness_p = abs(skewness_p)
plot(p, skewness_p, type="l", ylab="Skewness(abs)")
title("p의 변화에 따른 왜도(절대값)의 변화")
p[which(skewness_p == min(skewness_p))]
#원자료
stem(price)
#재표현자료
stem(price^3.4)
par(mfrow=c(1,2))
boxplot(price, col="papayawhip", ylim=c(1000000, 2000000), ylab="price(won)")
title("Boxplot of Samsung Electronics' stock price(2016)")
boxplot(price^3.4, col="papayawhip", ylab="price(won)")
title("Boxplot of Re-expressed Samsung Electronics' stock price(2016)")
s1 = abs(skewness(price))
s2 = abs(skewness(log(price)))
s3 = abs(skewness(price^0.5))
s4 = abs(skewness(price^-1))
s5 = abs(skewness(price^-0.5))
s6 = abs(skewness(price^2))
s = c(s1, s2, s3, s4, s5, s6)
s = round(s, 2)
s = as.character(s)
par(mfrow=c(2,3))
boxplot(price, col="papayawhip", ylab="price(won)")
title("original")
text(x= 0.8, y= 1700000, labels= paste("왜도:",s[1]), cex=2)
boxplot(log(price), col="papayawhip", ylab="price(won)")
title("log")
text(x= 0.8, y= 14.35, labels= paste("왜도:",s[2]), cex=2)
boxplot(price^0.5, col="papayawhip", ylab="price(won)")
title("p=0.5")
text(x= 0.8, y= 1300, labels= paste("왜도:",s[3]), cex=2)
boxplot(-price^-1, col="papayawhip", ylab="price(won)")
title("p=-1")
text(x= 0.8, y= -6.0e-07, labels=paste("왜도:",s[4]), cex=2)
boxplot(-price^-0.5, col="papayawhip", ylab="price(won)")
title("p=-0.5")
text(x= 0.8, y= -0.000765, labels= paste("왜도:",s[5]), cex=2)
boxplot(price^2, col="papayawhip", ylab="price(won)")
text(x= 0.8, y= 2.8e+12, labels= paste("왜도:",s[6]), cex=2)
title("p=2")