Week 01: Parameter Learning

해당 수업에서는 첫 번째 Machine Learning 알고리즘이면서, 빈번하게 사용되는 gradient descent 알고리즘을 학습 한다.

Gradient Descent Algorithm

Cost Fuction이 다음과 같이 있다고 가정하다.

당연히 이러한 Cost Fucntion의 값이 최소가 되는 θ0와 θ1을 우리는 찾고 싶을 것이다.

최적화된 Cost Function을 찾는 방법을 Gradient Descent 알고리즘이라 한다.

• Start with some θ0,θ1
• Keep Changing θ0,θ1 to reduce J(θ0,θ1) untile we hopfully end up at a minimum

아래와같이 최소 지역을 찾을 수고 있고,

아래와 같이 약간 옆으로 시작할경우 local minimum 문제에 빠질 수도 있다.

수식으로 Gradient descent algorithm을 표현하면 아래와 같다.
구현상의 중요점은 θ0,θ1두 개가 모두 한번에 업데이트 되어야 한다는 것이다.
그리고 temp를 이용해서 저장하고 그것을 반영해야한다.

Gradient Descent Algorithm Intuition

좀 더 직관적으로 쉽게 이해해보자.
일단 단하나의 변수 θ1만 존재하는 Cost Function을 생각해보자. 아래와 같이 간단히 구할 수 있다.

추가로 learning rate에 따라서 최적화가 될 수도 있고 안될 수도 있다.

아래는 간단한 Quize이다.

이렇게 local minimum에서 빠져 나오지 못하는것은
자체적으로 점점 더 값이 작아지기 때문이다.
알파 값을 고저이켜도 결국 minimum 값에 converge 하게 된다.

이러한 gradient descent algorithm을 이용해서 어떠한 cost function도 최적화 할 수 있다.
단 local minimum의 문제는 여전히 존재 한다.

Gradient Descent for Linear Regression

실제 cost function인 linear regression에서의 measn sqaured error 를 생각해보자.
당연히 θ0,θ1으로 두개의 변수를 가지고 있다.
각각에 대해서 partial derivative (편미분)하면 아래와 같다.

실제 partialderivative (편미분)을 통해서 optimal 함수를 구하는 것은 알아냈다.
이제는 좀 더 실제적으로 local fucntion을 구하는것을 알아보자.

bow shape을 따른다 Linear regression의 cost function의 경우에는

아래와같은 모양을 convex fucntion이라고 하고 이것을 bow-shaped 이라고 인포멀리 부른다.
이러한 convex fucntion은 local minimum이 존재하지 않고, 오로지 global minimum만이 존재한다.

따라서 linear regression의 cost function은 convex이기 때문에 local minimum문제 없이 언제나 global minimum에 도달하게 된다.

좀 더 포멀하게 cost function을 최적화 하는 것을 생각해보자.
θ0=900,θ1=0.1의 값을 생각해보자.
이때 hypothesis function h(x) = -0.1x+900

그리고 이것을 반복하면 최종적으로 global minimum에 도달하게 된다.

그리고 이렇게 전체 Traning set을 가지고 gradient descent하는 것을 Batch Gradient Descent라 한다.

Quize

data(iris); library(ggplot2); library(caret)

## Loading required package: lattice

names(iris)

## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 
## [5] "Species"

table(iris$Species)

## 
##     setosa versicolor  virginica 
##         50         50         50

inTrain <- createDataPartition(y=iris$Species,
                              p=0.7, list=FALSE)
training <- iris[inTrain,]
testing <- iris[-inTrain,]
dim(training); dim(testing)

## [1] 105   5

## [1] 45  5

modlda = train(Species ~ .,data=training,method="lda")

## Loading required package: MASS

modnb = train(Species ~ ., data=training,method="nb")

plda = predict(modlda,testing); pnb = predict(modnb,testing)
table(plda,pnb)

##             pnb
## plda         setosa versicolor virginica
##   setosa         15          0         0
##   versicolor      0         17         0
##   virginica       0          1        12

Comparsion of results
see that just one value appears between the two classes appears to be not classified in the same way by the same way two algorithms but overall they perform very similarly.

equalPredictions = (plda==pnb)
qplot(Petal.Width,Sepal.Width,colour=equalPredictions,data=testing)

Boosting is one of the bset box classifier

Basic Idea
Take lots of (poosibly) weak predictors

Weight them and add them up

Get a stronger predictor

Start with a set of classifiers h1,…,hk

Boosting in R
Boosting can be used with any subset of classifiers One large subclass is gradient boosting R has multiple boosting libraries. Differences include the choice of basic classification functions and combination rules. gbm - boosting with trees. mboost - model based boosting ada - statistical boosting based on additive logistic regression gamBoost for boosting generalized additive models Most of these are available in the caret package

Wage example

library(ISLR); data(Wage); library(ggplot2); library(caret);

## Loading required package: lattice

Wage <- subset(Wage,select=-c(logwage))
inTrain <- createDataPartition(y=Wage$wage,
                              p=0.7, list=FALSE)
training <- Wage[inTrain,]; testing <- Wage[-inTrain,]

#Fit the model
modFit <- train(wage ~., method="gbm",data=training,verbose=FALSE)

## Loading required package: gbm
## Loading required package: survival
## 
## Attaching package: 'survival'
## 
## The following object is masked from 'package:caret':
## 
##     cluster
## 
## Loading required package: splines
## Loading required package: parallel
## Loaded gbm 2.1.1
## Loading required package: plyr

print(modFit)

## Stochastic Gradient Boosting 
## 
## 2102 samples
##   10 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 2102, 2102, 2102, 2102, 2102, 2102, ... 
## Resampling results across tuning parameters:
## 
##   interaction.depth  n.trees  RMSE      Rsquared   RMSE SD   Rsquared SD
##   1                   50      35.15759  0.3038776  1.530055  0.02487676 
##   1                  100      34.62096  0.3137957  1.439743  0.02232634 
##   1                  150      34.55009  0.3156775  1.394560  0.02166566 
##   2                   50      34.60094  0.3156275  1.455682  0.02485196 
##   2                  100      34.47184  0.3182034  1.374485  0.02300061 
##   2                  150      34.54167  0.3162264  1.387161  0.02253147 
##   3                   50      34.47797  0.3187507  1.440898  0.02512358 
##   3                  100      34.59459  0.3142979  1.374172  0.02372285 
##   3                  150      34.82958  0.3071383  1.363458  0.02276088 
## 
## Tuning parameter 'shrinkage' was held constant at a value of 0.1
## 
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## RMSE was used to select the optimal model using  the smallest value.
## The final values used for the model were n.trees = 100,
##  interaction.depth = 2, shrinkage = 0.1 and n.minobsinnode = 10.

qplot(predict(modFit,testing),wage,data=testing)

bootstrap aggregating의 short name이다.

어떤 복잡한 모델을 생성할때 여러번 겹쳐서 평균을 낸다면

bias와 variance를 동시에 캡쳐할 수 있다.

기본 아이디어

리샘플해서 모델을 만들고 리프리덱션해서 정확도를 계산해서

그것의 평균을 구하게 된다.

리샘플의 방법은 bootstrapping을 따르는 복원추출을 이용한다.

특징

비슷한 bias를 계속해서 누적하며

하지만 variance를 줄이는 효과가 있다.

이것은 non-linear function을 쓸때 가장 유용하다.

예제: Ozone data

install.packages("ElemStatLearn")
library(ElemStatLearn); data(ozone,package="ElemStatLearn")
ozone <- ozone[order(ozone$ozone),]
head(ozone)
    ozone radiation temperature wind
17      1         8          59  9.7
19      4        25          61  9.7
14      6        78          57 18.4
45      7        48          80 14.3
106     7        49          69 10.3
7       8        19          61 20.1

오존으로 정렬한것이다.

하려고하는것은 temperature를 ozone function으로 예측할려고 하는 것이다.

이제 수작업으로 bagging을 수행해보자.

10 by 155 짜리 matrix를 생성 한다.

10번동안 resampleling을 해서 데이터셋을 만들고

그것으로 smooth curve를 loess 함수를 이용해서 만들고

이것으로 prediction하고

평균을 구하게 된다.

ll <- matrix(NA,nrow=10,ncol=155)
for(i in 1:10){ #10번을 반복 수행함
    ss <- sample(1:dim(ozone)[1],replace=T) #복원 추출을 수행 한다.
    ozone0 <- ozone[ss,]; ozone0 <- ozone0[order(ozone0$ozone),] #데이터를 추출하고 ozone으로 정렬 한다
    loess0 <- loess(temperature ~ ozone,data=ozone0,span=0.2) # loess를 이용해서 예측 모델을 생성 한다.
    ll[i,] <- predict(loess0,newdata=data.frame(ozone=1:155)) 
    # 예측 모델을 이용해서 1~155사이의 오존을 예측한다. 
    #실제 오존은 168번까지 있다. 이것은 순차적이지 않으며, 중간 중가 결측값(NA)도 존재한다.
}

그리고 실제 Ozone과 temperature에 대한 것을 검은색 점으로 표현하고

예측한 값들을 선으로 연결해서 10개를 그리고

그 10개의 예측값을 편균으로해서 빨간색 선으로 그리면 아래와 같은 그래프를 얻을 수 있다.

plot(ozone$ozone,ozone$temperature,pch=19,cex=0.5)
for(i in 1:10){lines(1:155,ll[i,],col="grey",lwd=2)}
lines(1:155,apply(ll,2,mean),col="red",lwd=2)

위와 같이 bias는 높지만 변동은 더 줄어든 것을 볼 수 있다.

정리

bagging은 nonlinear model들에서 많이 사용된다.

tree에서 많이 사용된다. 이것을 확장한 것이 random forests 이다.

기억할것은 bagging이란 결국

데이터를 복원주출로 샘플링 한다음 비선형모델로 생성을 하고

그렇게 생성한 여러개의 모델을 평균화해서 좀더 variance(분산)이적은 부드러운 모델을 생성하는 것이다.

이러한 부드러운 모델은 각각 데이터 셋으로 모델링 했을 때는 나오기 어려운 것들이다. 즉 상대적으로 변동이 적은 그러나 bias는 높은 모델을 생성 하는 것이다. 그것이 bagging 이다.

Ensemble method: Bagging (bootstrap aggregating)

일반적인 ensemble 방법의 하나가 bootstrap aggregating인 Bagging 방법이다.

Leo Breiman에 의해서 1994년에 만들어 졌다.

bagging은 오리지널 트레이닝 셋으로 부터 bootstrap 샘플링 방법을 이용해서 여러개의 트레이닝 데이터를 생성하는 방법이다.

이렇게 생성된 데이터 셋은 하나의 모델을 생성하기 위해서 사용 된다.

상대적으로 Ensemble 기법들 중에서는 Bagging은 단순한 방법에 속하지만 꽤 잘 동작 한다.

특히 Decision Tree 알고리즘의 경우 약간의 입력 데이터의 변화에도 급격히 다른 결과값을 이끌어 내므로

Bagging으로 이것을 보완 할 수 있다.

ipred package는 bagged decision tree를 구현한 것이다.

resubstitution error는 한글로 재치환오류라고 하며 트레이닝 셋에서 훈련된 모델에서의 에러를 말한다.

다른 용어로는 In sample error라고도 한다.

당연히 수치는 긍정적으로 평가된다. 왜냐하면 같은 데이터로 훈련하고 테스트 하기 때문이다.

Bootstrap sampling

큰 데이터 셋을 추정하기 위해서 임의로 여러번 샘플링을 하는것을 말한다.

When this principle is applied to machine learning model performance, it implies the creation of several randomly-selected training and test datasets, which are then used to estimate performance statistics.

The results from the various random datasets are then averaged to obtain a final estimate of future performance.

어떤 점이 그럼 k-fold Cross Validation이랑 다른 것일까?

이 방법은 생성한 각각의 구획들이 단 한번만 나타나게 된다.

하지만 bootstrap의 방법은 복원 추출을 하기 때문에 여러번 나올 수 있다.

증명에 의하면 초기 데이터 셋의 수인 N이 엄청 크다면, bootstrap으로는 63.2%를 커버한다는 증명이 있다.

결국 36.8%는 트레이닝 셋이다.

original data의 size가 N 레코드라고 하자. 이때 bootstrap으로 샘플링을 수행하면 N번 수행하면 당연히 N개의 레코드가 있고

이러한 N개의 레코드는 원래 original data size인 N개 중 63.2%를 포함한다고 한다.

N번의 bootstrap sampling에서 단일 record가 선택 되어질 확률은 아래와 같다.

$$1-(1-\frac{1}{N})^{N}$$

이것은 asymptotically approaches $1-e^{-1}=0.632$에 근접한다.

이렇게 bootstrap sampling으로 포함되지 않은 것들은 test set으로 자동으로 포함된다.

이에반에 10-fold cv를 할경우 90%가 결국 트레이닝이 된다.

N이 충분히 크다면 bootstrap은 더 낮은 representative를 보인다.

당연히 63.2%의 트레이닝 셋을 가지고 테스트 하기 때문에 더 큰 트레이닝셋을 가지고 만든 모델보다 더 낮은 정확도를 가지게 된다.

최종적인 error rate의 계산식은 아래와 같다.

$$ error = 0.632 \times error_{test} + 0.368 \times error_{train}$$

cross validation 보다 더 좋은 bootstap의 장점은 더 작은 dataset에서 잘 동작하는 특성이 있다는 것이다.

또한 bootstrap sampling은 모델 정확도 향상의 효과가 있을 수도 있다

주어진 데이터로 부터 복원 표본을 구하는 작업을 여러 번 반복해 원하는 값을 추정 하는 것이다.

예를 들어 보면,

평균의 신뢰 구간을 구하는 경우를 생각해보자.

평균의 신뢰 구간은 평균이 속할 범위를 95% 신뢰도로 찾는 것이다.

부트스트래핑은 데이터로 부터 표본을 복원 추출로 구하고 이들의 평균을 기록하는 일을 반복한다.

이렇게 구한 평균들을 나열한 뒤 가장 작은 값으로 부터 2.5% 가장 큰 값으로 부터 97.5% 지점의 평균 두 개를 구한다.

그러면 두 평균 사이의 구간이 평균 95% 신뢰 구간이 된다.

실제 코드로 구현해 보자.

평균 30, 표준편차 3인 정규 분포로부터 1,000개의 난수를 만들고 이 난수로부터 역으로 평균을 추정해보자.

x <- rnorm(1000, mean=30, sd=3)
t.test(x)
	One Sample t-test
data:  x
t = 308, df = 999, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 29.72329 30.10447
sample estimates:
mean of x 
 29.91388 

t검정 결과 평균의 신뢰 구간이 (29.78891, 30.16317)로 구해졌다. 

library(foreach)
bootmean <- foreach(i=1:10000, .combine=c) %do% {
    return(mean(x[sample(1:NROW(x), replace=TRUE)]))
}

이제는 부트스트래핑으로 복원 추출을 통해 샘플을 구하고 평균을 구하는 작업을 10,000회 반복해보자.

평균의 신뢰구간은 이 평균들의 하위 2.5%, 상위 97.5%에 속하는 값의 구간이다.

bootmean <- sort(bootmean)
bootmean[c(10000*0.025, 10000*0.975)]
[1] 29.72578 30.10531

t검정과 유사한 값이 부트스트랩핑으로도 구해진것을 알 수 있다.

부트스트래핑을 사용한 모델의 정확도 평가도 마찬가지로 수행된다.

데이터의 일부를 트레이닝 데이터로, 일부를 테스트 데이터로 복원 추출을 통해서 구한다.

그리고 이를 수회 반복하여 정확도를 추정한다. 이것이 trainControl()에서 method="boot"가 하는 일이다.

.

이론적인 내용은 이전 포스트를 참조 한다.

Tree 방식의 Classification의 장단점은 아래와 같다.

Pros:

쉬운 이해

비선형 데이터에 대한 좋은 성능

Cons:

pruning이나 cross-validation을 하지 않으면 over-fitting의 문제가 발생함.

linear regression model같은 것들보다 불확실한 데이터 즉 트레이닝에서 발견되지 않았던 많이 다른 데이터를

추정하기아 어렵다.

결과가가 항상 특정 tree를 구성하는 node 값에의해서 결정되므로 트리에 구조에따라서 결과가 자주 변한다.

예제는 아래와 같다.

2010년 미국 대통령 선거당시 힐러리 클린턴과 버락 오바에 대한 것이다.

각각의 질의에 따라 어떻게 투표 결과가 이뤄 질지에대한 예측을 나타내는 트리이다.

예제: Iris Data

해당 코드는 caret package를 이용하며 해당 package의 party, rpart 방식을 이용해서 tree를 구성한 것이다.

다른 package로는 tree package가 존재 한다.

데이터를 불러온다.

data(iris); library(ggplot2)
names(iris)
[1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"

각 테이블의 양을 파악한다

table(iris$Species)
    setosa versicolor  virginica 
        50         50         50 

트레이닝 데이터와 테스트 데이터를 분할 한다.

inTrain <- createDataPartition(y=iris$Species,
                               p=0.7, list=FALSE)
training <- iris[inTrain,]
testing <- iris[-inTrain,]
dim(training); dim(testing)
[1] 251   5
[1] 33  5

그래프를 그려서 세개의 분포를 확인해 보자.

qplot(Petal.Width,Sepal.Width,colour=Species,data=training)

아래와 같이 세개의 종들이 나름 극명하게 나뉘는 것을 알 수 있다. 

이제 caret package를 이용해서 모델을 학습해서 수행해 보자.

~,의 의미는 모든 variables를 species를 prediction하기위해서 사용하겠다는 것이다.

library(caret)
modFit <- train(Species ~ .,method="rpart",data=training)
print(modFit$finalModel)
n= 117 
node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 117 76 virginica (0.34188034 0.30769231 0.35042735)  
  2) Petal.Length< 2.45 40  0 setosa (1.00000000 0.00000000 0.00000000) *
  3) Petal.Length>=2.45 77 36 virginica (0.00000000 0.46753247 0.53246753)  
    6) Petal.Width< 1.75 40  5 versicolor (0.00000000 0.87500000 0.12500000) *
    7) Petal.Width>=1.75 37  1 virginica (0.00000000 0.02702703 0.97297297) *

모델을 그래프로 표현한다. 좀 더 명확히 보기위해서 rattle package를 이용해서 그린다.

아래의 그림은 일반 plot()이용해서 그린 그래프 이다.

library(rattle)
fancyRpartPlot(modFit$finalModel)
#draw plot using nomal plot
plot(modFit$finalModel, uniform=TRUE, 
     main="Classification Tree")
text(modFit$finalModel, use.n=TRUE, all=TRUE, cex=.8)

fancyRpartPlot으로 그리면 아래와 같다.

normal plot으로 그리면 아래와 같다.

생성된 tree 모델을 이용해서 예측을 수행 한다.

predict(modFit,newdata=testing)
 [1] setosa     setosa     setosa     setosa     setosa    
 [6] setosa     setosa     setosa     setosa     setosa    
[11] versicolor versicolor versicolor versicolor versicolor
[16] versicolor versicolor versicolor versicolor versicolor
[21] versicolor versicolor versicolor versicolor virginica 
[26] virginica  virginica  virginica  virginica  virginica 
[31] virginica  virginica  virginica 
Levels: setosa versicolor virginica

Dimensionality Reduction을 수행 할경우 아래와 같은 장점이 존재 한다.

데이터 압축의 효과 발생 -> 낮은 메모리 사용 디스크 공간 활용 증대 -> 속도 향상

따라서 중요하 전처리 과정중 하나이다.

Dimensonlaity Reduction의 실체

두 개의 차원의 데이터가 cm와 inches로 나뉘어서 있다고 하자. 

높은 중복성을 가지는 데이터이다.

이것을 파란색 선으로 그어보면 새로운 일차원 데이터 구조로 만들 수 있음을 알수 있다.

2D를 1D로 바꾼 것이 된다.

또 다른 예제는 아래와 같다.

헬기 조종사의 적성을 평가할때 skill과 enjoying 두개의 요소가 있다고 하자. 두개의 관계도 대칭적이므로 일차원으로 표현 할 수 있다.

결국 아래와 같은 데이터 압축의 효과가 나타난다.

3D to 2D의 경우도 아래와 같다.

데이터 표현의 경우 아래와 같다. 50개의 차원을 2개의 차원으로 줄인 것이다.

D50 -> D2를 함으로써 데이터를 그래프로 표현 할 수 있게 되었다.

Principal Component Analysis (PCA)

블루 라인이 projection error 이다.

이것은 sum of squares of these little blue line segments is minimized로 구해진다.

즉 최소제곱법 합으로 구해진다.

PCA는 이러한 최소화된 project 방법을 찾는 알고리즘이다.

mean normalization at feature scaling 

마젠타 라인에서의 블루라인은 길기때문에 PCA를 이용해서 최적의 선을 찾아내는 것이 중요하다.

formulation으로 표현해보자.

2차원에서 1차원으로 줄일때의 필요한 방향 벡터 $u^{(1)}$를 찾는것이 중요하다.

아래와 같이 해당 벡터는 -일수도 있고 +일 수도 있다. 단위 벡터를 찾기만 하면 되지 방향은 중요하지 않다.

n차원을 k 차원으로 줄일려면 k개의 벡터를 찾는것이 된다. $u^{(1)}, u^{(2)},...,u^{(k)}$와 같이이다.

3차원 2차원 예제로 본다면 아래와 같이 k-vector들에 걸쳐 있는 linear subspace를 찾을 수 있게 된다.

PCA는 뭔가 Linear regression을 의미하는것 같지만 그것과 동일 하지는 않다.

겉으로는 비슷해 보일 지언정 정말로 다른 알고리즘이다.

선형회귀분석의 경우 블루라인의 제곱합을 최대한 작게 하는것이다.

즉 sum of squared errors (SSE)가 가장 작아야 한다.

블루라인은 Vertical 라인이다.

하지만 PCA의 경우 vertical이 아니다.

아래의 그래프와 같이 포인트 X와 read line에 대해서 orthogonal (직교)한 최단거리를 의미하게 된다.

linear regression에서는 y를 predict 하기 위해서 x를 사용하는 것이 된다.

PCA에서는 distinguished y, special variable y 따위는 존재하지 않는다. 

즉 $x_1, x_2, ... , x_n$까지 모두 같은 동일한 것으로 다루게 된다.

아래의 예제처럼 $x_1, x_2, x_3$의 세개의 features가 존재하는 분포가 있다고 하자.

3차원 분포를 2차원으로 줄이기 위해서는 $u^{(1)}, u^{(2)}$ 두개의 방향 벡터를 찾아야 한다.

이때 $x_1, x_2, x_3$ 세개의 feature들은 모두 같은 수준이지 특별히 y가 없다.

아래의 문제를 풀면

$u^{(1)}$는 $[\frac{-1/\sqrt{2}}{1/\sqrt{2}}]$의 값을 가진다.

주성분 분석은 데이터에 많은 변수가 있을 때 변수의 수를 줄이는 차원 감소(Dimensionality Reduction)기법 중 하나다.

PCA는 변수들을 주성분(Principal Component)이라 부르는 선형적인 상관관계가 없는 다른 변수들로 재표현한다.

무슨말인가 하면,

선형적인 관계가 있는 두 변수가 아래와 같이 있다고 가정해보자.

A = 2c+3

B= 3c

이다.

이때 Y = f(A,B)라는 어떤 모델이 있다고 하자.

이럴때 실제로 A,B 둘 사이에는 C라는 변수에 의한 선형적 관계가 있는데 둘 다 이용해서 Y를 예측하는 모델을 생성할 필요가 있을까?

이와 같을 때는 f(A,B)로 Y를 예측하기 보다는 그냥 f(C)로 Y를 예측하는것이 더 좋다.

Principal component들은 원 데이터의 분산(퍼짐 정도)을 최대한 보존하는 방법으로 구한다.

좀 더 자세한 내용은 아래의 이전 포스트를 참조하자.

practical machine learning (johns hopkins)

machine learning (standford)

R을 이용해서 실제로 적용해 보기 [R프로그래밍 실습서]

pincomp() 함수 또는 preComp()를 이용한다.

pincomp: 주성분 분석을 수행 한다.

아래와 같이 1:10을 저장한 x

약간의 noise를 추가한 y. 즉 x + noise = y 이다.

마지막으로 z는 x + y 에 noise를 추가한 것이다.

data frame으로 내용을 출력해보면 아래와 같다.

위 데이터들은 주성분들이 원 데이터의 분산 중 얼마만큼을 설명해 주는지를 알 수 있다.

Proportion of Variance 행을 보면 첫 번째 주성분은 데이터의 분산 중 92.86%를 설명해 주며,

두 번째 주성분은 데이터의 분산 중 7.1%를 설명함을 알 수 있다.

세 번째 주성분은 가장작은 0.03%의 분산을 설명 한다.

마지막 행의 Cumulative Proportion은 Proportion of Variance의 누적 값이다.

결국 원 데이터의 분산은 첫 번째와 두 번째 주성분에 의해 99.97%가 포함됨을 알 수 있다.

이들 두 주성분상의 좌표는 scores를 보고 구하면 된다.

결국 x,y,z의 데이터 분포는 그냥 2개의 차원으로 축소가 된것이다. 아래와 같은 좌표값을 가지게 된다.

알츠하이머(Alzheimer)를 이용한 예제 [Practical Machine Learning HW-2]

데이터를 불러온다.

위 데이터 셋에서 "IL"로 시작하는 트레이닝 셋 variables(변수)들에 대해서만 PCA를 수행한다.

preProcess()를 이용하며, PCA 수행결과 80%의 variance(분산)를 capture할 수 있는 변수들의 갯수를 결정해 보자.

먼저 grep과 정규표현으로 해당 열 데이터를 수집한다.

아래와 같이 caret package의 preProcess()를 이용해서 pca를 수행한다. 임계값은 80%이다.

실행 결과 총 12개의 IL 변수중 7개만으로 원래 데이터 분포의 80%수준의 분산을 설명 할 수 있음을 알 수 있다.

이전과 같이 princomp를 수행하면 아래와 같다.

summary(princomp(training[, IL_str]))

PCA가 정확도에 미치는 영향을 glm 모델을 만들어서 알아보자 [Practical Machine Learning HW-2]

PCA를 적용하지 않은 glm 모델

PCA를 적용해서 수행한 결과이다.

적용전 0.65

적용후 0.72

결국 dimentionality reduction이 정확도 향상에도 기여한다는 것을 알 수 있다.

TensorFlow 설치법과 PyDev 설치법은 이전에 다뤘다.

원할한 개발을 위해서 PyDev에 TensorFlow를 연결하는것을 다룬다.

우선 필자는 TensorFlow를 Virtualenv에다가 설치 했다. VirtualEnv를 설치하는 이유는 아래와 같다.

PyDev의 VirtualEnv 설정 연결

PyDev와 TensorFlow 그리고 virtualenv가 모두 설치되었어야 한다.

Eclipse 메뉴에서

Preferences -> PyDev > Interpreters > Python 메뉴에서

New를 눌러서 Browse를 통해서 설정한 virtualenvs를 등록해준다.

등록할때 모든 종속 라이브러릴 같이 등록해준다. 팝업이 나오면 다 선택된 상태로 등록하면 된다.

프로젝트 생성 및 "Hellow Tensor"실행

아래와 같이 프로젝트를 생성 할 때 interpreter를 생성한 virtualenv용 python 환경으로 잘 선택해 줘야 한다.

그다음 Hello tensorflow 코드를 실행하면 정상적으로 실행 되는것을 알 수 있다.

'''
Created on Nov 17, 2015

@author: root
'''
import tensorflow as tf
hello = tf.constant('Hello, TensorFlow!')
sess = tf.Session()
print sess.run(hello)

a = tf.constant(10)
b = tf.constant(32)
print sess.run(a+b)

IntelliJ 기반의 Pycharm

설치방법: 이전포스트

실행 모습

아래와 같이 Tensorflow가 있는 virtualenv를 interpreter로 설정을 해주어야 한다.

웹기반 Jupyter로 실행하기

설치방법: 이전포스트

2015.11.9일에 공개된 오픈프로젝트 TensorFlow는 Google에서 공개한 DeepLearning을 위한 2세대 system이다.

TensorFlow의 내용은 이전 Post를 참조 한다.

이 시스템의 설치와 실행을 다룬다.

Python VirtualEnv에 기반한 설치 방법

필자가 사용하는 Desktop PC의 환경은 아래와 같다.

설치환경

필자는 VirtualEnv-based Installation을 수행한다(VirtualEnv 설치).

VirtualEnv 환경을 python 자체가 지원하므로 고립된 수행 환경을 생성할 수 있는 이점이 있다. 문제 발생시에도 모든 사람이 같은 상황이므로 해결방법을 그대로 적용할 수 있는 장점이 있다.

#Next, set up a new virtualenv environment. To set it up in the directory ~/tensorflow, run:
$ virtualenv --system-site-packages ~/tensorflow
$ cd ~/tensorflow

tensorflow는 이제 가상환경이 저장되는 디렉터리가 된다.

--system-site-packages 옵션의 의미는 가상환경에게 global site-package의 접근 권한을 허용 하는 것이다.

아래의 경로에 해당 가상환경에만 적용되는 site-package들이 저장되어 있다.

아래의 경로에는 해당 가상환경에서 사용 할 수 있는 python의 종류이다.

virtualenv를 활성화 시킨다.

$ source (설치경로)/bin/activate  # If using bash
source ./tensorflow/bin/activate #나의 경우 

비활성화는 deactivate라고 치면된다.

만약 csh 쉘을 쓴다면 아래의 명령어이다.

$ source bin/activate.csh  # If using csh

최종적으로 아래처럼 프롬프트 창의 이름 앞에 (tensorflow)가 들어가는 것을 볼 수 있다.

(tensorflow)root@jemin-virtual-machine:

실행된 virtualenv의 환경에서 TensorFlow를 설치한다.

설치 종류는 CPU-only version과 GPU-enabled version이 있으므로 자신이 원하는것을 선택하자.

필자는 일단 CUDA sdk를 설치하지 않았으므로 CPU-only를 설치했다.

# For CPU-only version
$ pip install https://storage.googleapis.com/tensorflow/linux/cpu/tensorflow-0.5.0-cp27-none-linux_x86_64.whl

# For GPU-enabled version (only install this version if you have the CUDA sdk installed)
$ pip install https://storage.googleapis.com/tensorflow/linux/gpu/tensorflow-0.5.0-cp27-none-linux_x86_64.whl

이렇게까지 하면 tensorFlow는 설치가 된것이다.

이제 간단한 예제를 실행해보자.

Hello, TensorFlow

Hello, TensorFlow
Python 2.7.6 (default, Jun 22 2015, 17:58:13) 
[GCC 4.8.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import tensorflow as tf
>>> hello = tf.constant("Hello, TensorFlow!")
>>> sess = tf.Session()
I tensorflow/core/common_runtime/local_device.cc:25] Local device intra op parallelism threads: 8
I tensorflow/core/common_runtime/local_session.cc:45] Local session inter op parallelism threads: 8
>>> print sess.run(hello)
Hello, TensorFlow!
>>> a = tf.constant(10)
>>> b = tf.constant(32)
>>> print sess.run(a+b)
42
(tensorflow)root@jemin-virtual-machine:~/tensorflow# 

# Make 100 phony data points in NumPy.
x_data = np.float32(np.random.rand(2, 100)) # Random input
y_data = np.dot([0.100, 0.200], x_data) + 0.300

# Construct a linear model.
b = tf.Variable(tf.zeros([1]))
W = tf.Variable(tf.random_uniform([1, 2], -1.0, 1.0))
y = tf.matmul(W, x_data) + b

# Minimize the squared errors.
loss = tf.reduce_mean(tf.square(y - y_data))
optimizer = tf.train.GradientDescentOptimizer(0.5)
train = optimizer.minimize(loss)

# For initializing the variables.
init = tf.initialize_all_variables()

# Launch the graph
sess = tf.Session()
sess.run(init)

# Fit the plane.
for step in xrange(0, 201):
    sess.run(train)
    if step % 20 == 0:
        print step, sess.run(W), sess.run(b)

# Learns best fit is W: [[0.100  0.200]], b: [0.300]

I tensorflow/core/common_runtime/local_device.cc:25] Local device intra op parallelism threads: 8
I tensorflow/core/common_runtime/local_session.cc:45] Local session inter op parallelism threads: 8
Hello, TensorFlow!
42
0 [[ 0.3138563   0.52770293]] [-0.0221094]
20 [[ 0.17014106  0.30057704]] [ 0.20599112]
40 [[ 0.12187628  0.23032692]] [ 0.27126268]
60 [[ 0.10676718  0.20918694]] [ 0.29121917]
80 [[ 0.10208294  0.2027912 ]] [ 0.29731771]
100 [[ 0.10063919  0.20084961]] [ 0.29918078]
120 [[ 0.10019579  0.20025893]] [ 0.29974979]
140 [[ 0.1000599   0.20007896]] [ 0.2999236]
160 [[ 0.10001831  0.2000241 ]] [ 0.29997668]
180 [[ 0.1000056   0.20000736]] [ 0.29999286]
200 [[ 0.1000017   0.20000222]] [ 0.29999784]

일단 gitHub에서 TensorFlow 소스코드를 복제(clone)해와야 한다.

tensorflow)root@jemin-virtual-machine: git clone https://github.com/tensorflow/tensorflow.git

정상적으로 다운이 완료 되면, 예제를 실행 해보자.

이 예제는 "LeNet -5 - like convolutional MNIST model"이라고 한다.

(tensorflow)$ cd tensorflow/models/image/mnist
(tensorflow)$ python convolutional.py

실행결과

Initialized!
Epoch 0.00
Minibatch loss: 12.054, learning rate: 0.010000
Minibatch error: 90.6%
Validation error: 84.6%
Epoch 0.12
Minibatch loss: 3.285, learning rate: 0.010000
Minibatch error: 6.2%
Validation error: 7.0%
Epoch 0.23
Minibatch loss: 3.473, learning rate: 0.010000
Minibatch error: 10.9%
Validation error: 3.7%
Epoch 0.35
Minibatch loss: 3.221, learning rate: 0.010000
Minibatch error: 4.7%
Validation error: 3.2%
Epoch 0.47
Minibatch loss: 3.193, learning rate: 0.010000
Minibatch error: 4.7%
Validation error: 2.7%
Epoch 0.58
Minibatch loss: 3.301, learning rate: 0.010000
Minibatch error: 9.4%
Validation error: 2.5%
Epoch 0.70
Minibatch loss: 3.203, learning rate: 0.010000
Minibatch error: 6.2%
Validation error: 2.8%
Epoch 0.81
Minibatch loss: 3.022, learning rate: 0.010000
Minibatch error: 4.7%
Validation error: 2.6%
Epoch 0.93
Minibatch loss: 3.131, learning rate: 0.010000
Minibatch error: 6.2%
Validation error: 2.1%
Epoch 1.05
Minibatch loss: 2.954, learning rate: 0.009500
Minibatch error: 3.1%
Validation error: 1.6%
Epoch 1.16
Minibatch loss: 2.854, learning rate: 0.009500
Minibatch error: 0.0%
Validation error: 1.8%
Epoch 1.28
Minibatch loss: 2.825, learning rate: 0.009500
Minibatch error: 1.6%
Validation error: 1.4%
Epoch 1.40
Minibatch loss: 2.938, learning rate: 0.009500
Minibatch error: 7.8%
Validation error: 1.5%
Epoch 1.51
Minibatch loss: 2.767, learning rate: 0.009500
Minibatch error: 0.0%
Validation error: 1.9%
Epoch 1.63
Minibatch loss: 2.771, learning rate: 0.009500
Minibatch error: 3.1%
Validation error: 1.4%
Epoch 1.75
Minibatch loss: 2.844, learning rate: 0.009500
Minibatch error: 4.7%
Validation error: 1.2%
Epoch 1.86
Minibatch loss: 2.694, learning rate: 0.009500
Minibatch error: 0.0%
Validation error: 1.3%
Epoch 1.98
Minibatch loss: 2.650, learning rate: 0.009500
Minibatch error: 0.0%
Validation error: 1.5%
Epoch 2.09
Minibatch loss: 2.667, learning rate: 0.009025
Minibatch error: 1.6%
Validation error: 1.4%
Epoch 2.21
Minibatch loss: 2.658, learning rate: 0.009025
Minibatch error: 1.6%
Validation error: 1.2%
Epoch 2.33
Minibatch loss: 2.640, learning rate: 0.009025
Minibatch error: 3.1%
Validation error: 1.2%
Epoch 2.44
Minibatch loss: 2.579, learning rate: 0.009025
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 2.56
Minibatch loss: 2.568, learning rate: 0.009025
Minibatch error: 0.0%
Validation error: 1.2%
Epoch 2.68
Minibatch loss: 2.554, learning rate: 0.009025
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 2.79
Minibatch loss: 2.503, learning rate: 0.009025
Minibatch error: 0.0%
Validation error: 1.2%
Epoch 2.91
Minibatch loss: 2.487, learning rate: 0.009025
Minibatch error: 0.0%
Validation error: 1.2%
Epoch 3.03
Minibatch loss: 2.463, learning rate: 0.008574
Minibatch error: 1.6%
Validation error: 1.2%
Epoch 3.14
Minibatch loss: 2.458, learning rate: 0.008574
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 3.26
Minibatch loss: 2.410, learning rate: 0.008574
Minibatch error: 0.0%
Validation error: 1.4%
Epoch 3.37
Minibatch loss: 2.496, learning rate: 0.008574
Minibatch error: 3.1%
Validation error: 1.3%
Epoch 3.49
Minibatch loss: 2.399, learning rate: 0.008574
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 3.61
Minibatch loss: 2.377, learning rate: 0.008574
Minibatch error: 0.0%
Validation error: 1.1%
Epoch 3.72
Minibatch loss: 2.333, learning rate: 0.008574
Minibatch error: 0.0%
Validation error: 1.1%
Epoch 3.84
Minibatch loss: 2.312, learning rate: 0.008574
Minibatch error: 0.0%
Validation error: 1.2%
Epoch 3.96
Minibatch loss: 2.300, learning rate: 0.008574
Minibatch error: 1.6%
Validation error: 1.1%
Epoch 4.07
Minibatch loss: 2.276, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 1.1%
Epoch 4.19
Minibatch loss: 2.250, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 4.31
Minibatch loss: 2.233, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 4.42
Minibatch loss: 2.217, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 4.54
Minibatch loss: 2.324, learning rate: 0.008145
Minibatch error: 3.1%
Validation error: 1.0%
Epoch 4.65
Minibatch loss: 2.212, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 4.77
Minibatch loss: 2.174, learning rate: 0.008145
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 4.89
Minibatch loss: 2.211, learning rate: 0.008145
Minibatch error: 1.6%
Validation error: 1.0%
Epoch 5.00
Minibatch loss: 2.193, learning rate: 0.007738
Minibatch error: 1.6%
Validation error: 1.0%
Epoch 5.12
Minibatch loss: 2.148, learning rate: 0.007738
Minibatch error: 3.1%
Validation error: 1.0%
Epoch 5.24
Minibatch loss: 2.153, learning rate: 0.007738
Minibatch error: 3.1%
Validation error: 1.0%
Epoch 5.35
Minibatch loss: 2.111, learning rate: 0.007738
Minibatch error: 1.6%
Validation error: 0.9%
Epoch 5.47
Minibatch loss: 2.084, learning rate: 0.007738
Minibatch error: 1.6%
Validation error: 0.8%
Epoch 5.59
Minibatch loss: 2.054, learning rate: 0.007738
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 5.70
Minibatch loss: 2.043, learning rate: 0.007738
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 5.82
Minibatch loss: 2.134, learning rate: 0.007738
Minibatch error: 3.1%
Validation error: 1.0%
Epoch 5.93
Minibatch loss: 2.006, learning rate: 0.007738
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 6.05
Minibatch loss: 2.048, learning rate: 0.007351
Minibatch error: 3.1%
Validation error: 0.9%
Epoch 6.17
Minibatch loss: 1.988, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 1.1%
Epoch 6.28
Minibatch loss: 1.957, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 6.40
Minibatch loss: 1.971, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 6.52
Minibatch loss: 1.927, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 6.63
Minibatch loss: 1.912, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 6.75
Minibatch loss: 1.901, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 6.87
Minibatch loss: 1.886, learning rate: 0.007351
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 6.98
Minibatch loss: 1.894, learning rate: 0.007351
Minibatch error: 1.6%
Validation error: 1.0%
Epoch 7.10
Minibatch loss: 1.859, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 7.21
Minibatch loss: 1.844, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 7.33
Minibatch loss: 1.836, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 1.0%
Epoch 7.45
Minibatch loss: 1.887, learning rate: 0.006983
Minibatch error: 3.1%
Validation error: 0.9%
Epoch 7.56
Minibatch loss: 1.808, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 7.68
Minibatch loss: 1.822, learning rate: 0.006983
Minibatch error: 1.6%
Validation error: 0.9%
Epoch 7.80
Minibatch loss: 1.782, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 7.91
Minibatch loss: 1.772, learning rate: 0.006983
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.03
Minibatch loss: 1.761, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.15
Minibatch loss: 1.773, learning rate: 0.006634
Minibatch error: 1.6%
Validation error: 0.9%
Epoch 8.26
Minibatch loss: 1.742, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.38
Minibatch loss: 1.744, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.49
Minibatch loss: 1.719, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 8.61
Minibatch loss: 1.700, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 8.73
Minibatch loss: 1.700, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 8.84
Minibatch loss: 1.801, learning rate: 0.006634
Minibatch error: 1.6%
Validation error: 0.8%
Epoch 8.96
Minibatch loss: 1.666, learning rate: 0.006634
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 9.08
Minibatch loss: 1.666, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 9.19
Minibatch loss: 1.649, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 9.31
Minibatch loss: 1.676, learning rate: 0.006302
Minibatch error: 1.6%
Validation error: 0.8%
Epoch 9.43
Minibatch loss: 1.626, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 9.54
Minibatch loss: 1.621, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 9.66
Minibatch loss: 1.606, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.8%
Epoch 9.77
Minibatch loss: 1.596, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.9%
Epoch 9.89
Minibatch loss: 1.602, learning rate: 0.006302
Minibatch error: 0.0%
Validation error: 0.9%
Test error: 0.8%

실행 시간은 아래와 같다.

- MNIST -

real 26m26.382s

user 77m58.714s

sys 19m21.539s

가상화 환경에서 CPU만 가지고 처리해서 오래걸린것 같다.

그래도 i7 4세대 프로세서, ram 24gb, ssd pro 256gb에서 실행한것인데, 대략 Android 5.0 컴파일 시간의 반 정도의 시간이 걸린것 같다.

python 설치방법 -기본

수동 설치 (비추천)

#download using the following command:
cd ~/Downloads/
wget https://www.python.org/ftp/python/2.7.9/Python-2.7.9.tgz

#Extract and go to the directory:
tar -xvf Python-2.7.9.tgz
cd Python-2.7.9

#Now, install using the command you just tried:
./configure
make
sudo make install

우분투에서 python3을 default로 잡지 않은 이유는 다음과 같다.

1) 우분투는 python 2.x의 스크립트를 많이 이용하고 있다.

2) python 2와 python3은 서로를 완전히 보완하지 못한다. 버전간의 충돌이 있다.

3) python 2.7은 지속적으로 bugfix가 되며, 2020년까지는 지원을 해준다고 한다.

4) 대부분의 python script는 명시적으로 #!/usr/bin.env python 이라고만 명시하는데, 만약 /usr/bin/python을 python3으로 변경 한다면, #/usr/bin/env python2라고 명시한 스크립트를 제외하고는 python2에서만 동작하는 스크립트가 문제가 생기게 된다.

python 설치방법 - Virtualenv 

# On Linux:
$ sudo apt-get install python-pip python-dev python-virtualenv
#활성화 방법
source (설치경로)/bin/activate #일반적인 경우
source ./tensorflow/bin/activate #나의 경우
#비활성화 방법
deactivate

IDE 설정 - PyDev (Eclipse)

필자는 PyDev를 사용한다.

하지만 다른 IDE가 궁금하다면 10 Best Python IDE for Developers를 읽어보자.

필자가 PyDev를 선택한 이유는 간단하다. Eclipse 기반의 Plugin으로써 Eclipse의 개발 경험을 그대로 활용할 수 있기 때문이다.

그리고 commercial IDE는 사용하고 싶지 않기 때문이다.

설치방법

Help -> Install New Software -> Add

Name: Pydev and Pydev Extensions

Location: http://pydev.org/updates

• Caret package
• Data slicing
• Training options
• Plotting predictions
• Basic preprocessing
• Covariate creation
• Preprocessing with principal components analysis
• Predicting with Regression
• Predicting with Regression Multiple Covariates

Predicting with Regression Multiple Covariates

다루는 내용은 두가지로,

첫 번째는 mutiple covariates를 이용해서 regression model을 생성하는 방법을 다루며,

두 번째는 mutiple covariates들 중에서 무엇을 선택해야하는지를 알아내는 방법이다.

Example: Predicting Wage

# wage data
library(ISLR); library(ggplot2); library(caret);
data(Wage); Wage <- subset(Wage,select=-c(logwage))
summary(Wage)

# get training / test set
inTrain <- createDataPartition(y=Wage$wage,
                               p=0.7, list=FALSE)
training <- Wage[inTrain,]; testing <- Wage[-inTrain,]
dim(training); dim(testing)

# feature plot to know how those are related to each other
featurePlot(x=training[,c("age","education","jobclass")],
            y = training$wage,
            plot="pairs")

# plot age versus wage
qplot(age,wage,data=training)

# to figure out the reason to make outlier
# color point by jobclass, most of point are blue.
qplot(age,wage,colour=jobclass,data=training)

# color education
# also explaine, advanced degree is outlier
qplot(age,wage,colour=education,data=training)

# fit a linear model
modFit<- train(wage ~ age + jobclass + education,
               method = "lm",data=training)
finMod <- modFit$finalModel
print(modFit)

plot(finMod,1,pch=19,cex=0.5,col="#00000010")
qplot(finMod$fitted,finMod$residuals,colour=race,data=training)
plot(finMod$residuals,pch=19)

# predicted versus truth in test set
pred <- predict(modFit, testing)
qplot(wage,pred,colour=year,data=testing)

# all covariates
modFitAll<- train(wage ~ .,data=training,method="lm")
pred <- predict(modFitAll, testing)
qplot(wage,pred,data=testing)

• Caret package
• Data slicing
• Training options
• Plotting predictions
• Basic preprocessing
• Covariate creation
• Preprocessing with principal components analysis
• Predicting with Regression
• Predicting with Regression Multiple Covariates

Predicting with Regression

선형적이면 구현하기 쉽고 이해가 쉽다.

하지만 비선형이면 정확도가 너무 떨어진다.

Example: old faithful eruptions

# ppredicting with regression
library(caret);data(faithful); set.seed(333)
inTrain <- createDataPartition(y=faithful$waiting,
                               p=0.5, list=FALSE)
trainFaith <- faithful[inTrain,]; testFaith <- faithful[-inTrain,]
head(trainFaith)

  eruptions waiting
1     3.600      79
3     3.333      74
5     4.533      85
6     2.883      55
7     4.700      88
8     3.600      85

# showing graph
plot(trainFaith$waiting,trainFaith$eruptions,pch=19,col="blue",xlab="Waiting",ylab="Duration")

데이터가 선형이므로 Linear Regression을 적용하면 좋다.

# making model
lm1 <- lm(eruptions ~ waiting,data=trainFaith)
summary(lm1)
Call:
lm(formula = eruptions ~ waiting, data = trainFaith)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.26990 -0.34789  0.03979  0.36589  1.05020 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.792739   0.227869  -7.867 1.04e-12 ***
waiting      0.073901   0.003148  23.474  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.495 on 135 degrees of freedom
Multiple R-squared:  0.8032,	Adjusted R-squared:  0.8018 
F-statistic:   551 on 1 and 135 DF,  p-value: < 2.2e-16

# model fit
plot(trainFaith$waiting,trainFaith$eruptions,pch=19,col="blue",xlab="Waiting",ylab="Duration")
lines(trainFaith$waiting,lm1$fitted,lwd=3)

만들어진 모델을 확인한다.

그리고 값을 prediction을 해보면 아래와 같다.

# predict a new value
coef(lm1)[1] + coef(lm1)[2]*80
(Intercept) 
   4.119307 
newdata <- data.frame(waiting=80)
predict(lm1,newdata)
       1 
4.119307 

테스트와 트레이이닝에 대해서 각각 그려보면 아래와 같다.

# plot predictions - training and test
par(mfrow=c(1,2))
plot(trainFaith$waiting,trainFaith$eruptions,pch=19,col="blue",xlab="Waiting",ylab="Duration")
lines(trainFaith$waiting,predict(lm1),lwd=3)
plot(testFaith$waiting,testFaith$eruptions,pch=19,col="blue",xlab="Waiting",ylab="Duration")
lines(testFaith$waiting,predict(lm1,newdata=testFaith),lwd=3)

각각의 모델의 에러를 계산해보면 아래와 같다. 계산 방식은 RMSE를 이용하며 트레이닝과 테스트에 대해서 각각 수행 했다.

> # Calculate RMSE on training
> sqrt(sum((lm1$fitted-trainFaith$eruptions)^2))
[1] 5.75186
> # Calculate RMSE on test
> sqrt(sum((predict(lm1,newdata=testFaith)-testFaith$eruptions)^2))
[1] 5.838559

이러한 구현을 caret을 이용해서 할경우 매우 쉽게 할 수 있다.

> #same process with caret
> modFit <- train(eruptions ~ waiting,data=trainFaith,method="lm")
> summary(modFit$finalModel)

Call:
lm(formula = .outcome ~ ., data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.26990 -0.34789  0.03979  0.36589  1.05020 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.792739   0.227869  -7.867 1.04e-12 ***
waiting      0.073901   0.003148  23.474  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.495 on 135 degrees of freedom
Multiple R-squared:  0.8032,	Adjusted R-squared:  0.8018 
F-statistic:   551 on 1 and 135 DF,  p-value: < 2.2e-16

• Caret package
• Data slicing
• Training options
• Plotting predictions
• Basic preprocessing
• Covariate creation
• Preprocessing with principal components analysis
• Predicting with Regression
• Predicting with Regression Multiple Covariates

Covariate creation

covariates are sometimes called predictors and features

There are two levels of covariate creation, or feature creation

Level 1: From raw data to covariate

raw data takes the form of an image, text file, or a website.

That kind of information is very hard to build a predictive model around when you haven't summarized the information in some useful way in to either a quantitative or qualitative variable.

e-mail 내용이 아래와 같이 있다고 하자.

그냥 저 이메일을 처리할 수는 없기 때문에 2개의 feature를 끌어 낸다.

captialAve는 당연히 모두 대문자이므로 100% = 1 이다.

particular word frequency를 고려하는데, you가 있을 수 있다. 2이다.

number of dollar 의 경우 8개이다.

Level 2: Transforming tidy covariates

sqaure를 한다.

Level 1, Raw data -> covariates

Depends heavily on application

The balancing act is summarization vs. information loss

Examples:

Text files: frequency of words, frequency of phrases (Google ngrams), frequency of capital letters.

Images: Edges, corners, blobs, ridges (computer vision feature detection))

Webpages: Number and type of images, position of elements, colors, videos (A/B Testing)

People: Height, weight, hair color, sex, country of origin.

The more knowledge of the system you have the better the job you will do.

When in doubt, err on the side of more features

Can be automated, but use caution!

Level 2, Tidy covariates -> new covariates

More necessary for some methods (regression, svms) than for others (classification trees).

Should be done only on the training set

The best approach is through exploratory analysis (plotting/tables)

New covariates should be added to data frames

dummy variable

removing zero covariates

Spline basis

# spline basis
#data
library(ISLR); library(caret); data(Wage);
inTrain <- createDataPartition(y=Wage$wage,
                               p=0.7, list=FALSE)
training <- Wage[inTrain,]; testing <- Wage[-inTrain,]

library(splines)
bsBasis <- bs(training$age,df=3) 
bsBasis
#linear model
lm1 <- lm(wage ~ bsBasis,data=training)
plot(training$age,training$wage,pch=19,cex=0.5)
points(training$age,predict(lm1,newdata=training),col="red",pch=19,cex=0.5)

Notes and further reading

# Correlated predictors
library(caret); library(kernlab); data(spam)
inTrain <- createDataPartition(y=spam$type,
                               p=0.75, list=FALSE)
training <- spam[inTrain,]
testing <- spam[-inTrain,]

M <- abs(cor(training[,-58]))
diag(M) <- 0
which(M > 0.8,arr.ind=T)
       row col
num415  34  32
direct  40  32
num857  32  34
direct  40  34
num857  32  40
num415  34  40

# drawing graph
names(spam)[c(34,32)]
plot(spam[,34],spam[,32])

# we could rotate the plot
X <- 0.71*training$num415 + 0.71*training$num857
Y <- 0.71*training$num415 - 0.71*training$num857
plot(X,Y)

#PCA
smallSpam <- spam[,c(34,32)]
prComp <- prcomp(smallSpam)
plot(prComp$x[,1],prComp$x[,2])

• Caret package
• Data slicing
• Training options
• Plotting predictions
• Basic preprocessing
• Covariate creation
• Preprocessing with principal components analysis
• Predicting with Regression
• Predicting with Regression Multiple Covariates

Basic pre-processing 

어떤 머신러닝 알고리즘을 적용하던지간에 anomaly value를 제거하는 것은 중요하다.

그래프를 그리던지해서 noise 제거하는것은 정확히 해당 데이터로 부터 right signal을 추출해서 모델을 생성하도록 돕기

때문에 정확도 항샹에 큰 영향을 미치게 된다.

아래와 같은 spam 데이터 예제에대해서 일부만 빈도가 높을 수도 있다. 이러면 문제가 발생한다.

평균과 표준편차를 계산해 보면 아래와 같다.

1.1 standardizing - manual

$$E(X)=m,\quad \sigma (X)=\sigma,\quad Z=\frac { X-m }{ \sigma  }$$라두면, 평균과 표준편차의 성질을 이용하여 아래와 같이 전개 된다.

$$E(Z)=E\left( \frac { X-m }{ \sigma  }  \right) =E=\left( \frac { 1 }{ \sigma  } X-\frac { m }{ \sigma  }  \right) =\frac { 1 }{ \sigma  } E(X)-\frac { m }{ \sigma  } =0$$ 

$$\sigma(Z)=\sigma\left( \frac { X-m }{ \sigma  }  \right) =\sigma=\left( \frac { 1 }{ \sigma  } X-\frac { m }{ \sigma  }  \right) =\frac { 1 }{ \sigma  } \sigma(X)=1$$

이것을 training set에 대해서 표준화(standardizing)을 수행하면 아래와 같은 결과를 얻을 수 있다.

> #Standardizing
> trainCapAve <- training$capitalAve
> trainCapAveS <- (trainCapAve  - mean(trainCapAve))/sd(trainCapAve) 
> mean(trainCapAveS)
[1] -5.008709e-18
> sd(trainCapAveS)
[1] 1

test set에 standardizing을 적용할 때에 test set에 대한 mean과 sd를 적용하는 것이 아닌,

training set에서의 mean과 sd 값을 적용 해야 한다. 

> #standardizing - test set
> testCapAve <- testing$capitalAve
> testCapAveS <- (testCapAve  - mean(trainCapAve))/sd(trainCapAve) 
> mean(testCapAveS)
[1] 0.04085023
> sd(testCapAveS)
[1] 1.131577

당연히 mean과 sd는 0과 1이 아니다. 하지만 그것에 근접한 값이 나오기를 바랄 뿐이다.

1.2 standardizing - using preProcess in caret package

Pre-processing transformation (centering, scaling etc.) can be estimated from the training data and applied to any data set with the same variables.

preObj <- preProcess(training[,-58],method=c("center","scale"))
trainCapAveS <- predict(preObj,training[,-58])$capitalAve
mean(trainCapAveS)
[1] -5.008709e-18
> sd(trainCapAveS)
[1] 1

여기서 58번째 컬럼을 제거하는 이유는 해당 컬럼만 numeric이 아닌 factor type 이기 때문이다.

method=c("center","scale")의 의미는 아래와 같다.

이전과 같이 거의 zero에 가까운 mean과 1의 sd를 얻을 수 있다.

test set의 경우 당연히 training set의 preObj을 사용해서 처리 해야한다.

> testCapAveS <- predict(preObj,testing[,-58])$capitalAve
> mean(testCapAveS)
[1] 0.04085023
> sd(testCapAveS)
[1] 1.131577

1.3 standardizing - preProcess argument

> set.seed(32343)
> modelFit <- train(type ~.,data=training,
+                   preProcess=c("center","scale"),method="glm")
There were 27 warnings (use warnings() to see them)
> modelFit
Generalized Linear Model 

3451 samples
  57 predictor
   2 classes: 'nonspam', 'spam' 

Pre-processing: centered (57), scaled (57) 
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 3451, 3451, 3451, 3451, 3451, 3451, ... 
Resampling results

  Accuracy   Kappa     Accuracy SD  Kappa SD 
  0.9201069  0.832212  0.01353308   0.0260575

2. Standardizing - Box-Cox transforms

새로운 standardizing 방법이 있다.

center, scale 방법은 익스트림 값을 제거하는 스타일의 방법이다.

하지만 여기서 소개하는 box-cox 방식은 normal data 처럼 만들어주는 standardizing 방법이다.

# Standardizing box cox transforms
preObj <- preProcess(training[,-58],method=c("BoxCox"))
trainCapAveS <- predict(preObj,training[,-58])$capitalAve
par(mfrow=c(1,2)); hist(trainCapAveS); qqnorm(trainCapAveS)

오른쪽은 이전에 standardizing을 수행하지 않고 histogram을 그린것이다. 보는것과 같이 특정 구간에만 빈도가 몰린것을 알 수 있다. 하지만 box-cox를 적용한다음에는 분포가 고른것을 볼 수 있다. 

Q-Q plot을 통해서 해당 standardizing이 얼마만큼 normal distribution과 일치 하는지를 알 수 있다.

변위치(Quantiles) 0에 많은 데이터가 분포하므로 완벽히 45도 기울기의 일직선이 아느므로 완벽하지는 않다. 하지만 거의 유사하다. 

3. Standardizing - Imputing data

missing data를 처리하기 위한 전처리 방법이다. NA(missing data), 결측치가 있을경우 대부분의 통계 함수들이나 prediction algorithm들은 동작 자체를 하지 않는다. 따라서 반드시 처리해줘야하는 것들 중 하나이다.

이때 이러한 missing data를 impute(덮어씌우기)하는 알고리즘들을 배워보자.

k-nearest neighbor's imputation이 쉬운 방법중 하나이다.

코드는 아래와 같다.

#standardizing - imputing data
set.seed(13343)

# Make some values NA
training$capAve <- training$capitalAve
selectNA <- rbinom(dim(training)[1],size=1,prob=0.05)==1
training$capAve[selectNA] <- NA

# Impute and standardize
preObj <- preProcess(training[,-58],method="knnImpute")
capAve <- predict(preObj,training[,-58])$capAve

# Standardize true values
capAveTruth <- training$capitalAve
capAveTruth <- (capAveTruth-mean(capAveTruth))/sd(capAveTruth)

결과는 아래와 같다.

> quantile(capAve - capAveTruth)
           0%           25%           50%           75%          100% 
-3.5036289018 -0.0005761963  0.0006479362  0.0011937179  0.1832545934 
> quantile((capAve - capAveTruth)[selectNA])
           0%           25%           50%           75%          100% 
-3.5036289018 -0.0175633552  0.0009851735  0.0159918619  0.1832545934 
> quantile((capAve - capAveTruth)[!selectNA])
           0%           25%           50%           75%          100% 
-0.8523414202 -0.0004922746  0.0006475486  0.0011696849  0.0016053800 

Notes and future reading

Training and test must be processed in the same way

Test transformations will likely be imperfect

Especially if the test/training sets collected at different times

Careful when transforming factor variables!

더 읽어볼 거리

preprocessing-caret: http://topepo.github.io/caret/preprocess.html

2015.11.9일 구글이 자시들의 Machine Learning Platform TensorFlow를 공개 했다.

공식 사이트: http://tensorflow.org/

GitHub: https://github.com/tensorflow/tensorflow

ResearchBlog: http://googleresearch.blogspot.kr/

이것은 구글이 2011년부터 내부에서 사용해오던 DeepLearning Infrastructure인 "DistBelief"의 후속 버전이다.

DistBelief는 지금까지 음성인식, 지메일 스팸 방지, 유튜브 동영상 썸네일 제작, 유트부 영상에서 고양이 인식 등의 딥러닝을 적용한 구글의 서비스와 연구에 사용된 것이다.

이번에 공개된 TensorFlow는 이러한 DistBelief의 성능과 사용편의성을 개선한 2세대 머신 러닝 시스템이라고 할 수 있다.

GoogleResearch - TensorFlow, Google's latest machine learning system, open sourced for everyone.

작성자: Jeff Dean (Senior Google Fellow) and Rajat Monga (Technical Lead)

2011년 부터 사용된 DistBelief 구글의 DeepLearning 플랫폼은 Find "cat", speech recognition, Google Photos, Large Scale Visual Recognition Challegne in 2014, and automated image captioning 그리고 DeepDread 에서 성공적인 결과를 보였다.

하지만, DistBelief는 너무나 Google의 내부 구조와 결집되어 있기 때문에 그동안 공개하기가 어려웠다.

2015.11.9일 구글은 마침내 2세대 머신러닝 시스템인 TensorFlow를 공개한다. 이것은 그동안의 단점을 모두 보완 했으며, flexible, portable, easy to use, and completely open source이다.

또한 구글의 Technical report에 따르면 TensorFlow는 특정 BanchMark에 대해서 2배 가량 DistBreief보다 빠르다고 한다.

이제 non-learning tool과 learning too로 양분화 대며, 모든 non-learning tool들은 점차 대체될 것이다.

바야흐로 새로운 시대가 열렸다.

TensorFlow의 특징

일단은 구글이 사용하는 버전과 다르게 단일 컴퓨터에서 다중 GPU를 지원하여 동작하는것 같다.

Jeff Dean의 말에 따르면 이것은 기존의 Deep Learning library인 Torch, Caffe, Theano의 장점을 가져온 것이라고 한다.

즉, Torch와 Theano 처럼 research idea를 구현하기 쉽고 Caffe처럼 실제 product에 구현하기 쉽게 만든것이 Tensorflow이라는 것이다.

참고문헌: http://www.wired.com/2015/11/google-open-sources-its-artificial-intelligence-engine/

• Caret package
• Data slicing
• Training options
• Plotting predictions
• Basic peprocessing
• Covariate creation
• Preprocessing with principal components analysis
• Predicting with Regression
• Predicting with Regression Multiple Covariates

• RMSE = Root mean squared error
• RSquared = $R^2$ from regression models

• Accuracy = Fraction correct
• Kappa = A measure of concordance

library(ISLR); library(ggplot2); library(caret);
data(Wage)
summary(Wage)

      year           age               sex                    maritl    
 Min.   :2003   Min.   :18.00   1. Male  :3000   1. Never Married: 648  
 1st Qu.:2004   1st Qu.:33.75   2. Female:   0   2. Married      :2074  
 Median :2006   Median :42.00                    3. Widowed      :  19  
 Mean   :2006   Mean   :42.41                    4. Divorced     : 204  
 3rd Qu.:2008   3rd Qu.:51.00                    5. Separated    :  55  
 Max.   :2009   Max.   :80.00                                           
                                                                        
       race                   education                     region    
 1. White:2480   1. < HS Grad      :268   2. Middle Atlantic   :3000  
 2. Black: 293   2. HS Grad        :971   1. New England       :   0  
 3. Asian: 190   3. Some College   :650   3. East North Central:   0  
 4. Other:  37   4. College Grad   :685   4. West North Central:   0  
                 5. Advanced Degree:426   5. South Atlantic    :   0  
                                          6. East South Central:   0  
                                          (Other)              :   0  
           jobclass               health      health_ins      logwage     
 1. Industrial :1544   1. <=Good     : 858   1. Yes:2083   Min.   :3.000  
 2. Information:1456   2. >=Very Good:2142   2. No : 917   1st Qu.:4.447  
                                                           Median :4.653  
                                                           Mean   :4.654  
                                                           3rd Qu.:4.857  
                                                           Max.   :5.763  
                                                                          
      wage       
 Min.   : 20.09  
 1st Qu.: 85.38  
 Median :104.92  
 Mean   :111.70  
 3rd Qu.:128.68  
 Max.   :318.34  

위의 것을 그래프로 그려보자.

featurePlot()은 caret package에 있는 것이다.

Wrapper for Lattice Plotting of Predictor Variables

#drawing Feautre plot
featurePlot(x = training[,c("age","education","jobclass")],
            y = training$wage,
            plot="pairs")

재미있는 그래프 이다. wage에 대해서 즉 y값에 대해서 age, education, jobclass와의 관계를 한눈에 볼 수 있다.

다소 복잡해 보이지만 찬찬히 보면 한눈에 들어오는것을 알 수 있다.

• Caret package
• Data slicing
• Training options
• Plotting predictions
• Basic peprocessing
• Covariate creation
• Preprocessing with principal components analysis
• Predicting with Regression
• Predicting with Regression Multiple Covariates

caret package

해당 페키지는 R에서 Machine Learning을 수행하도록 도와주는 많은 package들 중의 하나이다.

하지만 여타 package들과 다르게 간편하고 강력하다는 장점이 있다.

■ caret의 기능들

1) 프리프로세싱 (데이터 클린)

2) 데이터 분할

createDataPartition

createResample

createTimeSlices

3) 트레이닝과 테스트 함수

train

predict

4) 생성한 모델들간의 비교

confusionMatrix

caret에서 지원하는 Machine Learning Algorithms

• Linear discriminant analysis 

• Regression

• Naive Bayes

• Support vector machines

• Classificiation and Regression tree

• Random forests
• Boosting
• etc.

왜 caret 인가? 수많은 machine learning 구현 package들 중에서

caret이 설정해서 사용하기 제일 편하다.

아래의 표처럼 각각의 페키지들은 서로다른 알고리즘을 사용할 때 마다 설정을 재각각 다르게 해줘야 한다.

다시보는 SPAM 예제

Data Splitting 

Fit a model

train 함수를 이용해서 machine learning을 수행한다.

prediction할 것은 type이고

tilde (~)와 dot(.) 은 모든 data frame에 있는 데이터들은 prediction을 위해서 사용된다는 것을 의미한다.

그리고 어떤 데이터셋을 트레이닝을 위해서 사용할지 알려준다.

data=training

마지막으로 어떤 method를 사용할지를 알려준다.

method="glm"

생성된 모델 정보는 아래와 같다.

모델의 계수(coeeficients)들을 확인하면 아래와 같다.

모델을 이용해서 prediction할 수 있다.

predict() 함수를 이용하는데, 이때 통일된 인자를 넘겨 주게된다. 이러한 통일성이 caret의 장점이다.

인자로는 생성된 model과 prediction에 사용될 데이터 set을 넘겨 준다.

모델의 평가는 confusion matrix을 이용해서 한다.

많은 statistical information들이 자동으로 계산 되는것을 볼 수 있다.

추가적인 참고 자료들

공식 사이트: http://topepo.github.io/caret/index.html

튜토리얼 자료1: http://www.edii.uclm.es/~useR-2013/Tutorials/kuhn/user_caret_2up.pdf

튜토리얼 자료2: https://cran.r-project.org/web/packages/caret/vignettes/caret.pdf

caret 논문 (journal of statistical software): file:///C:/Users/justin/Downloads/v28i05.pdf

Data Slicing

data slicing에 대한 이야기 이다.

뭘하던 간에 데이터 슬라이싱 부터 수행 해야 한다.

caret에서 지원하는 Data Splitting function들은 아래와 같다.

각각에 대해서 알아보자.

■ createDataPartition

아래의 createDataPartion() 함수를 이용해서 데이터를 분할 하게 된다.

y 변수에 분할할 데이터를 입력하고

p에서 몇 %를 분할 할지를 설정한다.

그럼 최종적으로 inTrain에는 index 값이 저장 되어 진다.

■ createFolds

cross validation방법중 k-folds 방법으로 data와 training set을 분할 하는 방법을 소개한다.

y는 결과값을 저장하고 있는 vector 이다.

k는 folds할 숫자를 나타낸다.

list = TRUE로 설정할 경우 결과를 list로 저장 하게 된다.

returnTrain = TRUE로 설정할 경우 folds의 결과를 training을 위해서 사용할지 여부를 결정 한다.

생성된 10-folds의 데이터 길이를 각각 확인하면 아래와 같다.

각각의 fold의 데이터 길이는 유사한것을 알 수 있다.

각 fold가 포함하는 데이터를 10개씩만 확인해보면 아래와 같다.

test set을 생성하는 방법은 아래와 같다.

■ createResample

완전히 full cross validation하는 것이 아니라 resampling 또는 bootstrapping을 하고 싶을 때가 있다.

설정방법은 아래와 같다.

y는 이전과 같이 결과 데이터를 저장하고 있는 vector 이다.

time은 얼마만큼 resampling을 할지를 저장해 준다.

list, vector, matrix 원하는 type을 설정해 준다.

결과를 확인해보면 k-folds 방법과 다르게 중복 데이터가 존재하는 것을 알 수 있다.

■ createTimeSlice

시계열 데이터를 slicing 할 필요가 있을 수도 있다.

기상예보를 예측하려고 할때를 말한다.

아래와 같이 time series 데이터의 경우 

window 크기를 20으로 설정한다음 그것으로 다음 10개를 예측 하는 식으로 동작하게 된다.

이렇게 순차적으로 데이터를 분할하게 된다.

즉, (트레이닝 20개, 테스트 10)이 조합이 반복적으로 수행 된다.

이론적인 부분은 아래의 이전 포스트를 확인한다.

Confusion Matrix

가장 쉽게 생성할 수 있는 방법은 table()을 이용 하는 것이다.

vec1 <- factor(c("ham","ham","spam","ham"))
vec2 <- factor(c("spam","ham","spam","ham"))
> table(vec1,vec2)
      vec2
vec1   ham spam
  ham    2    1
  spam   0    1

좀더 상세한 결과를 원한다면, gmodels package의 CrossTable() function을 사용하자.

default로써 CrossTable() output으로써 proportion값을 포함하고 있다.

> class(sms_test_pred)
[1] "factor"
> class(sms_raw_test$type)
[1] "factor"
> CrossTable(sms_test_pred, sms_raw_test$type)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  1390 

 
              | sms_raw_test$type 
sms_test_pred |       ham |      spam | Row Total | 
--------------|-----------|-----------|-----------|
          ham |      1202 |        31 |      1233 | 
              |    16.109 |   106.250 |           | 
              |     0.975 |     0.025 |     0.887 | 
              |     0.996 |     0.169 |           | 
              |     0.865 |     0.022 |           | 
--------------|-----------|-----------|-----------|
         spam |         5 |       152 |       157 | 
              |   126.514 |   834.437 |           | 
              |     0.032 |     0.968 |     0.113 | 
              |     0.004 |     0.831 |           | 
              |     0.004 |     0.109 |           | 
--------------|-----------|-----------|-----------|
 Column Total |      1207 |       183 |      1390 | 
              |     0.868 |     0.132 |           | 
--------------|-----------|-----------|-----------|

위 데이터를 가지고 직접 precision, recall, accuracy, error rate 등을 계산 하면 된다.

Beyond accuracy, other measures of performance

자동 계산을 위해서는 또 다른 package인 caret을 사용 한다.

해당 package는 classification과 regression training에 대한 정보를 담고 있으며, Max Kuhn에 의해서 만들어 졌다.

preparing, training, evaluating, and visualizing machine learning models and data 등을 지원 한다.

Journal of statistical software에 "Building predictive models in R using the caret package."을 보면 상세한 내용을 알 수 있다.

아래와 같이 사용하며, 단지 positive label이 무엇인지만 정의해 주면 된다. 

> library(caret)
필요한 패키지를 로딩중입니다: lattice
필요한 패키지를 로딩중입니다: ggplot2
> confusionMatrix(sms_test_pred, sms_raw_test$type, positive = "spam")
Confusion Matrix and Statistics

          Reference
Prediction  ham spam
      ham  1202   31
      spam    5  152
                                          
               Accuracy : 0.9741          
                 95% CI : (0.9643, 0.9818)
    No Information Rate : 0.8683          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.8795          
 Mcnemar's Test P-Value : 3.091e-05       
                                          
            Sensitivity : 0.8306          
            Specificity : 0.9959          
         Pos Pred Value : 0.9682          
         Neg Pred Value : 0.9749          
             Prevalence : 0.1317          
         Detection Rate : 0.1094          
   Detection Prevalence : 0.1129          
      Balanced Accuracy : 0.9132          
                                          
       'Positive' Class : spam    

이제 일반적인 statistics를 알아보자.

■ The kappa statistic

■ Leave One Out Cross - Validation

cross validation 방법의 일종으로써

총 N개의 샘플을 생성 한다.

이때 N개는 해당 데이터 셋의 샘플들의 총 수를 말한다

즉, 100개의 데이터 샘플이 있다면 N은 100이다.

이에대해서 이제 N개의 모델을 생성 한다. 그리고 각각의 모델을 생성할때

하나의 샘플씩을 제외한다. 즉 N개의 모델은 N-1개의 샘플 수로 만들어지는 것이다.

그리고 각각의 모델에 대해서 그 제외한 샘플 1개를 이용해서 test를 수행한다.

그렇게 나온 test 결과 N개에대한 평균 값을 구하게 된다.

MSE가 측정 방법이라고 한다면 아래의 식과 같다.

$MSE_{ i }=(y_{ 2 }-\hat { y_{ 2 } } )^{ 2 }$

구하려고 하는 성능의 척도는 결국 아래와 같다.

$CV_{(n)}=\frac{1}{n}\sum_{i=1}^{n}{MSE_{i}}$

LOOCV의 장점은 결국 모든 샘플에 대해서 다 한번씩은 test하는 것이 되므로

어떠한 랜덤성도 존재하지 않게 된다. 따라서 굉장히 stable한 결과를 얻을 수 있다.

단점으로는 N개의 모델을 만들어서 N번 테스트 하게 되므로 computing time이 굉장히 오래 걸리게 된다.

또한 K-fold Cross validation에 비해서 bias의 비중이 높게 된다. 증 다양성을 포함하기 어렵게 된다.

아래는 작업을 수행한 코드이다.

# LOOCV 처리방법
library(ISLR)
glm.fit=glm(mpg~horsepower, data=Auto) #mpg를 horsepower에 대해서 linear regression을 힙니다.
library(boot) #cv.glm을 사용하기 위한 라이브러리입니다.
cv.err=cv.glm(Auto, glm.fit) #cv.glm함수를 통해 LOOCV를 시행합니다.
cv.err$delta #delta는 cross-validation 결과들을 담고 있습니다.

#linear model, polynomial model에 한해서 LOOCV에 간단한 형태의 공식을 사용할 수 있습니다.
loocv=function(fit){
    h=lm.influence(fit)$h
    mean((residuals(fit)/(1-h))^2)
}

loocv(glm.fit)

cv.error=rep(0,5)
degree=1:5 #1차부터 5차함수까지 fitting해 봅시다.
for(d in degree){
    glm.fit=glm(mpg~poly(horsepower,d), data=Auto) #poly함수를 이용해 d차 함수를 fit할 수 있습니다.
    cv.error[d]=loocv(glm.fit) #LOOCV 결과를 각 벡터에 저장합니다.
}
plot(degree,cv.error,type="b") #LOOCV 결과를 plot합니다.

5차원에서 가장 에러가 적은것을 알 수 있다.

자료출처: https://smlee729.github.io/r/machine%20learning/2015/03/19/1-loocv.html

ROC curves

Receiver Operating Characteristic Curves (ROC curves)

prediction algorithm의 좋은 정도를 측정하는 좋은 measure 이다.

Why a curve ?

In binary classification you are predicting one of two categories

Alive / Dead

Click on add / don't click

But your prediction are often quantitative

Probability of being alive

Prediction on a scale from 1 to 10

The cutoff you choose gives different results

같은 모델이라도 0.7 부터 alive로 설정한지 그 이하로 설정 할지에 따라서 결과가 달라지게 된다.

ROC curves

상세한내용 Wekipedia

사람들은 ROC curve를 이용하는데

X축은 1 - specificity 

다른 말로, P(FP)와 같다.

정리하면 아래와 같다.

$$1\quad -\quad specificity\quad =\quad 1\quad -\quad \frac { TN }{ (FP+TN) } \\ =>\quad \frac { FP }{ FP+TN } \quad =\quad Probability\quad of\quad FP\quad =\quad P(FP)$$

Y축은 sensitivity

다른말로 P(TP)와 같다.

정리하면 아래와 같다.

$$Sensitivity\quad =\quad \frac { TP }{ TP+FN } \quad =\quad P(TP)$$

이제 해당 그래프로 무엇을 하냐면

어떠한 curve를 하나 생성하고 curve를 구성하는 각각의 점들은 모두 하나의 cut off가 된다.

어떤 특별한 cut off에서 그 지점에서의 P(FP)와 P(TP)를 구할 수 있게 된다.

이러한 curve를 통해서 해당 Algorithm이 좋은 안 좋은지를 평가 할 수 있다.

즉 각각의 cutoff로 구성된 어떠한 curve를 생성해서 해당 알고리즘의 goodness를 평가 할 수 있다.

X축은 1 - Specificity 이므로 값이 커질 수록 specificity가 작아지는 것이다.

Y축은 Sensitivity 이므로 값이 클 수록 좋다.

x축 0.2는 specificity가 0.8이므로 매우 좋은것이다. 하지만 그것을 따라 위로 가보면 빨간색 알고리즘 (NetChop)에 대해서 0.4가 약간 넘는 Sensitivity를 가지는 것을 알 수 있다. 이것은 그다지 높은 수치가 아니다.

따라서 이 X와 Y축 둘 사이에는 trade-off 관계가 존재하는 것을 알 수 있다.

Area under the curve (AUC)

AUC = 0.5 라는 것은 random guessing과 다를바 없으므로 매우 안좋은 성능을 나타낸다.

0.5 보다 낮다는 것은 동전 던지기 (flip a coin) 보다 더 안좋다는 것이기 때문에 심각하다고 할 수 있다.

AUC = 1 라는 것은 perfect classifier

일반적으로 AUC가 0.8 이상이라면 좋다고 생각한다.

Cross validation

아래와 같은 목적으로 자주 사용되는 도구이다.

detecting relevant features for building model

estimating their parameters when doing machine Learning

Key idea

트레이닝 데이터 셋에대한 Accuracy에 대해서는 언제나 낙관적이다. 그다지 유효하지 않다.

왜냐하면 항상 best model이 선택 되어 지기 때문이다.

더 나은 평가 방법은 independent set (test set accuracy)를 이용해서 평가 하는 것이다.

하지만 test set을 언제나 사용 할 수 있는 것은 아니다. 즉 모델을 만들고 있을때나 데이터가 너무 작아서 test set을 설정하지 못할 때이다. test set은 반드시 단 한번만 적용 해야 한다.

결국 training set과 test set의 경계는 모호해지며 따라서 training set으로 independent set을 만들어서 평가하는 방법을

고안한 것이 cross validation 이다.

각각의 다른 방식의 cross validation을 알아보자.

Random sub-sampling

랜덤하게 데이터를 training과 testing으로 구분해서 수행 한다.

K-fold

K 번 나눠서 수행 한다.

그리고 각각의 결과를 평균화하여 최종적인 결과로 사용 한다.

Leave one out

하나를 빼고 하나를 넣고 이렇게 순차적으로 반복해서 수행 한다.

고려사항 (Considerations)

time series 데이터들은 "chunk"단위로 수행 해야 한다.

k-fold cross validation을 위해서는

k가 너무 크면 bias는 값이 작아진다. 하지만 variance는 커진다.

k가 너무 작다면 bias는 값이 커진다. 하지만 variance는 작아진다.

LOOCV는 좀 더 높은 variance를 K-fold CV에 비해서 가지게 된다. 

왜냐하면, Training set들이 LOOCV의 경우 좀더 overlap 되기 때문이다.

Random sampling must be done without replacement

Random sampling with replacement is the bootstrap

Underestimates of the error

Can be corrected, but it is complicated (0.632 Bootstrap)

If you cross-validate to pick predictors estimate you must estimate error on independent data.

• Prediction motivation
• What is prediction?
• Relative importance of steps
• In sample and out of sample errors
• Prediction study design
• Types of errors
• Receiver Operating Characteristic
• Cross validation
• What data should you use ?