R - 로지스틱 회귀분석

데이터 탐색

GRE, GPA, RANK이 입학(admission)에 어떤 영향을 주는지 로지스틱 회귀분석을 통해 분석한다.

library(aod)
library(ggplot2)

mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

# view the first few rows of the data
head(mydata)

# 데이터의 대략적인 분포 확인
summary(mydata)

# 데이터 구조 확인
str(mydata)

# 변수별 표준편차 확인
sapply(mydata, sd)

# contingency table : xtabs(~ x + y, data)
xtabs(~admit+rank, data=mydata)

> head(mydata)
  admit gre  gpa rank
1     0 380 3.61    3
2     1 660 3.67    3
3     1 800 4.00    1
4     1 640 3.19    4
5     0 520 2.93    4
6     1 760 3.00    2

> summary(mydata)
     admit             gre             gpa        rank   
 Min.   :0.0000   Min.   :220.0   Min.   :2.260   1: 61  
 1st Qu.:0.0000   1st Qu.:520.0   1st Qu.:3.130   2:151  
 Median :0.0000   Median :580.0   Median :3.395   3:121  
 Mean   :0.3175   Mean   :587.7   Mean   :3.390   4: 67  
 3rd Qu.:1.0000   3rd Qu.:660.0   3rd Qu.:3.670          
 Max.   :1.0000   Max.   :800.0   Max.   :4.000  

> str(mydata)
'data.frame':	400 obs. of  4 variables:
 $ admit: int  0 1 1 1 0 1 1 0 1 0 ...
 $ gre  : int  380 660 800 640 520 760 560 400 540 700 ...
 $ gpa  : num  3.61 3.67 4 3.19 2.93 3 2.98 3.08 3.39 3.92 ...
 $ rank : Factor w/ 4 levels "1","2","3","4": 3 3 1 4 4 2 1 2 3 2 ...

> sapply(mydata, sd)
      admit         gre         gpa        rank 
  0.4660867 115.5165364   0.3805668   0.9444602 

> xtabs(~admit+rank, data=mydata)
     rank
admit  1  2  3  4
    0 28 97 93 55
    1 33 54 28 12

GLM을 통한 로지스틱 회귀모형 구축

rank 변수를 factor 타입으로 변경시킨다.

glm 을통해 로지스틱회귀모형을 fitting시킨다. family='binomial' 인자를 통해, glm으로 로지스틱 회귀모형을 쓸 수 있다. (link function이 logit function이라는 의미)

mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")

summary(mylogit)

> summary(mylogit)

Call:
glm(formula = admit ~ gre + gpa + rank, family = "binomial", 
    data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6268  -0.8662  -0.6388   1.1490   2.0790  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.989979   1.139951  -3.500 0.000465 ***
gre          0.002264   0.001094   2.070 0.038465 *  
gpa          0.804038   0.331819   2.423 0.015388 *  
rank2       -0.675443   0.316490  -2.134 0.032829 *  
rank3       -1.340204   0.345306  -3.881 0.000104 ***
rank4       -1.551464   0.417832  -3.713 0.000205 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 458.52  on 394  degrees of freedom
AIC: 470.52

Number of Fisher Scoring iterations: 4

로지스틱 회귀모형 결과 해석

Call : 구축한 모형에 대해 다시 상기시켜준다.

Deviance Residuals : Deviance residual에 대한 정보를 알려주는데, model fitting이 잘 되었는지에 대한 measure이다. 이를 통해 모델이 잘 적합됐는지를 평가할 수 있다.

Coefficient : 회귀계수와 그것들의 표준편차, z-statistics(wals's z-statistics), p-value를 나타낸다. 위 결과에서는 모든 변수가 유의한 것을 알 수 있다. 로지스틱 회귀모형에서는 회귀계수가 변수가 한 단위 증가했을 때 log(odds)의 증가량으로 해석할 수 있다.

example

gre가 1증가할 때, admission의 log odds(non-admission에 대한)가 0.0022 증가한다.

gpu의 경우도 마찬가지로 해석한다.

더미변수인 rank의 경우 해석이 약간 다른데, 예를 들어 rank2의 회귀계수 -0.67은 rank1에서 rank2로 바뀌었을 때, log(odds)의 변화량이다. 즉, rank1에 비해 rank2가 admission에 안좋은 영향을 준다는 것을 알 수 있다.

회귀계수들의 신뢰구간 얻기

Log-likelihood를 통해 구하는 법

confint(mylogit)

> confint(mylogit)
Waiting for profiling to be done...
                    2.5 %       97.5 %
(Intercept) -6.2716202334 -1.792547080
gre          0.0001375921  0.004435874
gpa          0.1602959439  1.464142727
rank2       -1.3008888002 -0.056745722
rank3       -2.0276713127 -0.670372346
rank4       -2.4000265384 -0.753542605

주어진 회귀계수의 표준편차를 이용해 신뢰구간을 구하는 법

## CIs using standard errors
confint.default(mylogit)

> confint.default(mylogit)
                    2.5 %       97.5 %
(Intercept) -6.2242418514 -1.755716295
gre          0.0001202298  0.004408622
gpa          0.1536836760  1.454391423
rank2       -1.2957512650 -0.055134591
rank3       -2.0169920597 -0.663415773
rank4       -2.3703986294 -0.732528724

이 신뢰구간은 이렇게도 계산 가능하다.예를 들어, gre의 95% 신뢰구간은 아래와 같이 구한다.

> 0.002264-1.96*0.001094
[1] 0.00011976
> 0.002264+1.96*0.001094
[1] 0.00440824

Wald Test를 통해 범주형 변수의 overall effect 파악하기

rank 변수를 dummy 변수로 만들어, reference인 rank1과 비교해 어떤 effect가 있는지 분석하였다. 하지만 aod 패키지의 walt.test를 이용하면, rank의 overall effect를 파악할 수 있다. [b:회귀계수, Sigma:error term의 공분산행렬, Terms:rank변수가 있는 열] wald.test는 순서를 통해 해당 범주형 변수의 위치를 파악하기 때문에 순서를 잘 신경 써야 한다.

wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)

> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)
Wald test:
----------

Chi-squared test:
X2 = 20.9, df = 3, P(> X2) = 0.00011

카이스퀘어 테스트의 p-value가 0.00011이므로 rank 변수는 유의하다.

변수간의 회귀계수가 같은지 검정하기

wald test를 통해 rank가 admission에 유의한 효과가 있다는 것을 파악하였고, rank1이라는 reference에 비해 rank2와 rank2가 유의한 영향을 준다는 것도 확인하였다. 한가지 의문은 rank2, rank3, rank4 의 회귀계수가 동일한지에 대한 것이다. 회귀 계수가 동일하다면, rank1일 때에 비해 rank2, rank3, rank4의 효과가 동일하다는 것을 의미한다.

이 구문은 rank2의 효과와 rank3의 효과가 동일한지를 검정한다.

l <- cbind(0, 0, 0, 1, -1, 0)
wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)

> l <- cbind(0, 0, 0, 1, -1, 0)
> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)
Wald test:
----------

Chi-squared test:
X2 = 5.5, df = 1, P(> X2) = 0.019

검정 결과, rank2, rank3의 효과는 유의하다 다른 것으라 나타났다.

OR과 회귀계수의 관계

앞서 로지스틱 회귀분석에서의 회귀계수는 log odds의 증가량이라고 언급하였다. 그러면 회귀계수에 exponential을 취하면, OR의 증분이 된다.

exp(coef(mylogit))

> exp(coef(mylogit))
(Intercept)         gre         gpa       rank2       rank3       rank4 
  0.0185001   1.0022670   2.2345448   0.5089310   0.2617923   0.2119375 

해석 ex) gpa가 1증가하면 admission의 non-admission에 대한 OR이 2.23배 증가한다.

예측하기

1. 다른 변수들을 고정하고 rank가 변할 때 예측값의 변화를 보기

newdata1 <- with(mydata, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))

> newdata1
    gre    gpa rank
1 587.7 3.3899    1
2 587.7 3.3899    2
3 587.7 3.3899    3
4 587.7 3.3899    4

newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response") # probability 점추정값
newdata1

> newdata1
    gre    gpa rank     rankP
1 587.7 3.3899    1 0.5166016
2 587.7 3.3899    2 0.3522846
3 587.7 3.3899    3 0.2186120
4 587.7 3.3899    4 0.1846684

2. rank, gpa를 고정한 후, gre의 효과 보기

newdata2 <- with(mydata, data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100),
                                              4), gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))
newdata2

> newdata2
         gre    gpa rank
1   200.0000 3.3899    1
2   206.0606 3.3899    1
3   212.1212 3.3899    1
4   218.1818 3.3899    1
5   224.2424 3.3899    1
6   230.3030 3.3899    1
7   236.3636 3.3899    1
8   242.4242 3.3899    1
9   248.4848 3.3899    1
10  254.5455 3.3899    1

newdata3 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type = "link", 
                                    se = TRUE))
newdata3 <- within(newdata3, {
  PredictedProb <- plogis(fit)
  LL <- plogis(fit - (1.96 * se.fit))
  UL <- plogis(fit + (1.96 * se.fit))
})

> head(newdata3)
       gre    gpa rank        fit    se.fit residual.scale        UL        LL PredictedProb
1 200.0000 3.3899    1 -0.8114870 0.5147714              1 0.5492064 0.1393812     0.3075737
2 206.0606 3.3899    1 -0.7977632 0.5090986              1 0.5498513 0.1423880     0.3105042
3 212.1212 3.3899    1 -0.7840394 0.5034491              1 0.5505074 0.1454429     0.3134499
4 218.1818 3.3899    1 -0.7703156 0.4978239              1 0.5511750 0.1485460     0.3164108
5 224.2424 3.3899    1 -0.7565919 0.4922237              1 0.5518545 0.1516973     0.3193867
6 230.3030 3.3899    1 -0.7428681 0.4866494              1 0.5525464 0.1548966     0.3223773

이 때 type="link" 옵션을 주면, link scale로 예측값을 내준다. 이는 logit 변환을 하기전의 예측값이다. 이를 통해 예측한 확률의 신뢰구간을 구할 수 있다. 왜냐하면 이 link scale의 예측값 (fit 변수)의 표준편차를 쉽게 알 수 있기 때문이다. 따라서 link scale의 예측값의 신뢰구간을 구한후 이를 다시 logit 변환하여 예측 확률의 신뢰구간을 얻을 수 있다. 이를 구현한 것이 바로 위의 코드이다.

plogis 함수를 link scale의 예측값에서 확률 추정값 p를 구해준다. 즉 log(p/1-p) = Xb 의 방정식을 풀어서 p를 구한다. 예를 들어 첫번째 라인의 경우, log(p/(1-p)) = -0.811 의 방정식을 양변에 exponential을 취해 풀면, 대략 p=0.308이 나온다.

newdata3의 예측값 시각화하기

rank를 고정시켜놓고 gre가 admission 확률에 어떤 영향을 주는지를 시각화

ggplot(newdata3, aes(x = gre, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
                                                                    ymax = UL, fill = rank), alpha = 0.2) + geom_line(aes(colour = rank),
                                                                                                                      size = 1)

출처

https://stats.idre.ucla.edu/r/dae/logit-regression/

원래 데이터

해당 지역의 10000명당 기관의 수를 나타낸 데이터입니다. (출처-통계청)

           sgg_kor facility_type 2003fac 2004fac 2005fac 2006fac 2007fac
 4               계  상급종합병원  0.0086  0.0086 0.00852  0.0087  0.0087
6               계      종합병원  0.0494  0.0494 0.05054  0.0510  0.0530
8               계          병원  0.1645  0.1755 0.18450  0.1937  0.2127
10              계      요양병원  0.0139  0.0231 0.04120  0.0727  0.1200
14              계      치과병원  0.0211  0.0221 0.02517  0.0274  0.0311
18              계      한방병원  0.0311  0.0320 0.03024  0.0292  0.0288
22              계        조산원  0.0143  0.0129 0.01055  0.0103  0.0104
24              계        보건소  0.0473  0.0475 0.04750  0.0478  0.0481
26              계      보건지소  0.2601  0.2605 0.25859  0.2571  0.2602
28              계    보건진료소  0.3865  0.3881 0.38504  0.3843  0.3875
30              계    보건의료원  0.0035  0.0035 0.00345  0.0034  0.0035
36            서울  상급종합병원  0.0195  0.0195 0.01942  0.0193  0.0196
38            서울      종합병원  0.0418  0.0409 0.04079  0.0406  0.0412
40            서울          병원  0.1158  0.1168 0.12431  0.1304  0.1432
42            서울      요양병원  0.0029  0.0058 0.01068  0.0251  0.0461

Melt

Melt는 Column을 Row로 바꾸는 함수입니다. (Key=id에 지정, 나머지값들은 Row로 변환)

facility <- melt(facility, id=c("sgg_kor", 'facility_type'))

sgg_kor facility_type variable  value
 1      계  상급종합병원  2003fac 0.0086
 2      계      종합병원  2003fac 0.0494
 3      계          병원  2003fac 0.1645
 4      계      요양병원  2003fac 0.0139
 5      계      치과병원  2003fac 0.0211
 6      계      한방병원  2003fac 0.0311

reshape

Reshape함수를 통해 Row를 Column으로 바꿀 수 있습니다. (Key를 idvar로 놓고, 컬럼으로 바꾸고 싶은 변수를 timevar에 지정)

facility <- reshape(facility, idvar=c("sgg_kor", "variable"), timevar="facility_type", direction="wide")

  sgg_kor variable value.상급종합병원 value.종합병원 value.병원 value.요양병원 value.치과병원      
계         2003fac             0.0086         0.0494     0.1645         0.0139         0.02      
    서울  2003fac             0.0195         0.0418     0.1158         0.0029         0.0350         
  강남구  2003fac             0.0373         0.0560     0.2425             NA         0.1679         
  강동구  2003fac             0.0209         0.0417     0.1043             NA             NA         
  강서구  2003fac                 NA         0.0185     0.1668         0.0185         0.0185         
  관악구  2003fac                 NA         0.0190     0.0949             NA             NA             

R 랜덤포레스트 사용하기

머신러닝 분야에서 많이 쓰이는 예제 데이터셋 Pima Diabetes 자료를 통해 간단한 Random Forest를 사용하는 방법을 알아보겠습니다.

라이브러리 로드

library(MASS)
library(randomForest)
library(caret)

MASS : 예제 데이터셋이 있는 라이브러리

randomForest : randomForest 이용을 위한 라이브러리

caret : 성능 평가시 confusion matrix를 그리기 위한 라이브러리

데이터 불러오기

data("Pima.tr")
data("Pima.te")

R에서 데이터 불러오는 방법은 라이브러리가 로드된 상태에서 data("데이터이름")이라고 실행하면, 작업공간상에 데이터가 생성됩니다.

> head(Pima.tr)
  npreg glu bp skin  bmi   ped age type
1     5  86 68   28 30.2 0.364  24   No
2     7 195 70   33 25.1 0.163  55  Yes
3     5  77 82   41 35.8 0.156  35   No
4     0 165 76   43 47.9 0.259  26   No
5     0 107 60   25 26.4 0.133  23   No
6     5  97 76   27 35.6 0.378  52  Yes

RandomForest 모형 만들기

set.seed(10)
rf.fit = randomForest(type ~ npreg + glu + bp + skin + bmi + ped + age
                      , data=Pima.tr, mtry = floor(sqrt(7)), ntree = 500, importance = T)
rf.fit

randomForest 모형을 만드는 문법은 위와 같습니다. hyper parameter로 이 모형에서는 mtry와 ntree를 사용하였습니다. mtry는 각각의 tree마다 몇 개의 feature를 사용할 것인지를 정하는 것입니다. 골드스탠다다는 보통, regression의 경우 변수갯수/3, classification의 경우 sqrt(변수갯수)를 사용합니다. 이 데이터셋은 classification 문제이므로, sqrt를 사용하였습니다. ntree는 tree의 총 갯수를 의미합니다.

테스트 데이터 생성

test_x = Pima.te[c("npreg", "glu", "bp", "skin", "bmi", "ped", "age")]
test_y = Pima.te$type

test_y를 불러올 때는 $로 불러와야지 vector로 불러와지므로 아래 confusionMatrix가 잘 실행됩니다.

성능 평가

y_pred = predict(rf.fit, test_x)

confusionMatrix(y_pred, test_y)

Confusion Matrix and Statistics

          Reference
Prediction  No Yes
       No  192  46
       Yes  31  63
                                          
               Accuracy : 0.7681          
                 95% CI : (0.7189, 0.8124)
    No Information Rate : 0.6717          
    P-Value [Acc > NIR] : 7.678e-05       
                                          
                  Kappa : 0.455           
 Mcnemar's Test P-Value : 0.1106          
                                          
            Sensitivity : 0.8610          
            Specificity : 0.5780          
         Pos Pred Value : 0.8067          
         Neg Pred Value : 0.6702          
             Prevalence : 0.6717          
         Detection Rate : 0.5783          
   Detection Prevalence : 0.7169          
      Balanced Accuracy : 0.7195          
                                          
       'Positive' Class : No    

변수 중요도

> importance(rf.fit)
             No       Yes MeanDecreaseAccuracy MeanDecreaseGini
npreg  9.467923  1.762747             9.098049         9.308814
glu   14.543787 16.760373            20.263701        21.771324
bp     4.210273 -8.185459            -1.740227         7.920807
skin   1.370477  1.662921             2.029630         9.350372
bmi    3.431359  7.500978             7.200368        13.322145
ped    4.957824  4.666885             6.778956        13.479470
age   12.162707  6.341079            13.414806        14.528669

변수 중요도 플롯

17페이지에 random forest hyper-parameter 관련한 설명들이 나온다.

https://cran.r-project.org/web/packages/randomForest/randomForest.pdf

Random Forest 관련 튜토리얼

https://machinelearningmastery.com/tune-machine-learning-algorithms-in-r/

Keras 멀티 GPU 이용하기

keras 2.0.9 버전부터 멀티 GPU를 손쉽게 활용할 수 있게 업데이트되었습니다. 이번 포스팅은 기존 모델을 멀티 GPU를 사용하여 트레이닝하는 예제입니다.

1. keras 2.0.9 이상에서 지원하므로 2.0.9 이상으로 업데이트

pip install --upgrade keras

2. keras 모델 생성

일반적인 keras CNN을 예로 들어 설명하겠습니다. 대략 이런 cnn model이 있다고 합시다. (패키지 임포트문 생략)

def cnn_model():
    # create model
    model = Sequential()
    model.add(Convolution2D(16, 3, 3, border_mode='same', activation='relu', input_shape=(256,256,1)))
    model.add(BatchNormalization())
    model.add(MaxPooling2D(pool_size=(2, 2)))
    model.add(Convolution2D(32, 3, 3, border_mode='same', activation='relu'))
    model.add(BatchNormalization())
    model.add(MaxPooling2D(pool_size=(2, 2)))
    model.add(Convolution2D(32, 3, 3, border_mode='same', activation='relu'))
    model.add(BatchNormalization())
    model.add(MaxPooling2D(pool_size=(2, 2)))
    model.add(Convolution2D(64, 3, 3, border_mode='same', activation='relu'))
    model.add(MaxPooling2D(pool_size=(2, 2)))
    model.add(Convolution2D(64, 3, 3, border_mode='same', activation='relu'))
    model.add(BatchNormalization())
    model.add(MaxPooling2D(pool_size=(2, 2)))
    model.add(Flatten())
    model.add(Dense(64, activation='relu'))
    model.add(Dropout(0.5))
    model.add(Dense(64, activation='relu'))
    model.add(Dropout(0.5))
    model.add(Dense(3, activation='softmax'))

    return model

2.0.9 버전 이후부터 제공하는 multi_gpu_model을 임포트

from keras.utils.training_utils import multi_gpu_model

모델을 만들고 multi_gpu_model을 통해서 multi gpu를 활용한다는 선언을 합니다. gpus 파라미터에 사용하고 싶은 gpu 수를 설정합니다. 모델을 compile하고 fit하면 트레이닝이 됩니다.

model = cnn_model()
model = multi_gpu_model(model, gpus=4)
model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
model.fit(train_x, train_y, validation_data=(test_x, test_y), nb_epoch=20, batch_size=32, verbose=1)

트레이닝이 잘됩니다.

Train on 215 samples, validate on 107 samples
Epoch 1/20
215/215 [==============================] - 1s - loss: 2.0448 - acc: 0.4465 - val_loss: 1.0884 - val_acc: 0.2710
Epoch 2/20
215/215 [==============================] - 0s - loss: 1.5974 - acc: 0.4791 - val_loss: 1.1431 - val_acc: 0.1963
Epoch 3/20
215/215 [==============================] - 0s - loss: 1.2704 - acc: 0.5302 - val_loss: 1.2524 - val_acc: 0.4486
Epoch 4/20
215/215 [==============================] - 0s - loss: 1.1777 - acc: 0.4512 - val_loss: 1.3994 - val_acc: 0.5327
Epoch 5/20
215/215 [==============================] - 0s - loss: 1.0600 - acc: 0.5953 - val_loss: 1.4960 - val_acc: 0.5888
Epoch 6/20
215/215 [==============================] - 0s - loss: 1.0486 - acc: 0.6512 - val_loss: 1.5189 - val_acc: 0.5888
Epoch 7/20
215/215 [==============================] - 0s - loss: 0.9968 - acc: 0.6186 - val_loss: 1.5151 - val_acc: 0.5888
Epoch 8/20
215/215 [==============================] - 0s - loss: 1.0375 - acc: 0.6000 - val_loss: 1.4286 - val_acc: 0.5888
Epoch 9/20
215/215 [==============================] - 0s - loss: 0.9987 - acc: 0.6093 - val_loss: 1.3102 - val_acc: 0.5701
Epoch 10/20

서버에서 nvidia-smi를 쳐보면 이렇게 모든 gpu가 모두 활용되는 것을 볼 수 있습니다.

참고

https://www.pyimagesearch.com/2017/10/30/how-to-multi-gpu-training-with-keras-python-and-deep-learning/

우분투에 R과 Rstudio Server 설치

1. R 설치

sudo apt-get install r-base

2. R studio 설치

아래 링크에서 보고 따라하면 됨

https://www.rstudio.com/products/rstudio/download-server/

이건 1.1.383 version 기준

sudo apt-get install gdebi-core

wget https://download2.rstudio.org/rstudio-server-1.1.383-amd64.deb

sudo gdebi rstudio-server-1.1.383-amd64.deb

3. 웹브라우저 접속

http://<server-ip>:8787

위 처럼 접속해서 운영체제 아이디/패스워드를 통해 로그인

Python - sklearn LabelEncoder, OnehotEncoder 사용

python에서 범주형 변수를 인코딩 하기 위하여 더미변수를 만들거나, one hot encoding을 합니다. 선형회귀 같은 기본적인 통계 모형에서는 더미변수를 많이 쓰지만, 일반적인 머신러닝/딥러닝 접근법에서는 one hot encoding을 많이합니다. one hot encoding을 하기위해서는 직접 함수를 만들어 할 수도 있지만 패키지를 사용하면 편합니다. 가장 편하다고 생각하는 방법이 바로 sklearn 패키지를 통한 one hot encoding 인데요. 예를 들어 아래와 같은 info라는 이름의 pandas dataframe 이 있을 때, class2를 예측변수 y라고 생각하여 one hot encoding 하는 방법을 알아보겠습니다.

Code

from sklearn import preprocessing

label_encoder = preprocessing.LabelEncoder()
onehot_encoder = preprocessing.OneHotEncoder()

train_y = label_encoder.fit_transform(info['class2'])
train_y = train_y.reshape(len(train_y), 1)
train_y = onehot_encoder.fit_transform(train_y)

LabelEncoder 결과

[0 0 2 2 0 0 2 2 2 2 0 2 0 0 2 0 2 0 2 0 2 0 2 1 2 0 2 2 1 2 0 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 0 2 2 2 0 2 2 2 2 2 0 2 2 1 2 2 1 2 2 2 2 0 0 2 0 2 2 2 2 2

OnehotEncoder 결과 - 최종 train_y

  (0, 0)	1.0
  (1, 0)	1.0
  (2, 2)	1.0
  (3, 2)	1.0
  (4, 0)	1.0
  (5, 0)	1.0
  (6, 2)	1.0
  (7, 2)	1.0
  (8, 2)	1.0
  (9, 2)	1.0
  (10, 0)	1.0
  (11, 2)	1.0
  (12, 0)	1.0
  (13, 0)	1.0

몇 줄 되지 않는 코드로 이처럼 one hot encoding을 구현할 수 있습니다. one hot encoding을 하기 전에 label encoding을 하는 이유는 one hot encoder의 인풋으로 숫자형만 올 수 있기 때문입니다. label encoder의 결과로 문자형 변수가 숫자형 변수 범주형으로 변경되게 되고 이를 one hot encoder에 fit_transform 해주면, 이와 같이 one hot encoding된 결과를 얻을 수 있습니다. 이 때, train_y 변수는 sparse matrix가되어 프린트하면 위와같이 나타납니다. 이는 matrix의 인덱스와 그에 해당하는 value를 나타낸건데 매트릭스에 값을 대입해보면 labelencoder 결과에 onehot encoding이 적용된 것을 확인할 수 있습니다.

Keras- Tensorflow Backend에서 특정 디바이스 사용법

Keras에 tensorflow backend를 사용할 때, 기본적으로 gpu를 사용하여 돌아갑니다. gpu가 있는데 cpu를 사용할 필요는 없겠지만, gpu의 성능을 측정하거나 하는 목적을 위해 강제로 cpu를 사용하도록 하는 방법에 대해 알아보겠습니다.

keras + tensorflow에서 사용가능한 하드웨어 자원보기

아래와 같이 tensorflow의 device_lib이라는 함수를 통해 현재 사용가능한 하드웨어 디바이스들의 리스트를 볼 수 있습니다. 여기서 name을 통해 해당 하드웨어에 접근할 수 있습니다. 제가 사용하는 환경의 경우 1080ti 4개를 사용하는데, 이 디바이스들의 리스트가 나오는 것을 볼 수 있습니다.

from tensorflow.python.client import device_lib
print(device_lib.list_local_devices())

[name: "/cpu:0"
device_type: "CPU"
memory_limit: 268435456
locality {
}
incarnation: ....
, name: "/gpu:0"
device_type: "GPU"
memory_limit: 295043072
locality {
  bus_id: 1
}
incarnation: .....
physical_device_desc: "device: 0, name: GeForce GTX 1080 Ti, pci bus id: 0000:02:00.0"
, name: "/gpu:1"
device_type: "GPU"
memory_limit: 366346240
locality {
  bus_id: 1
}
...

특정 device 사용하여 모델 구축, 트레이닝하기

아래와 같이 keras.backend.tensorflow_backend를 K 라는 이름으로 불러오고 with 문을 통해 해당 환경에서 모델을 구축후, fitting을 하면 해당 디바이스 환경에서 모델을 트레이닝 할 수 있습니다. 주의할 점은 model = Sequential()을 통해 모델을 생성할 때부터 with문을 통해 해당 device를 사용한다는 것을 알려줘야 한다는 것입니다. 단순히 model.fit 만 with 아래에서 돌리면 with 문이 먹히지 않더군요..

import keras.backend.tensorflow_backend as K
with K.tf.device('/gpu:0'):
    model = Sequential()
    model.add(Dense(512, input_dim=28*28, activation='relu'))
    model.add(Dense(256, activation='relu'))
    model.add(Dense(128, activation='relu'))
    model.add(Dense(32, activation='relu'))
    model.add(Dense(10, activation='softmax'))

    model.compile(loss='categorical_crossentropy',
                  optimizer='rmsprop',
                  metrics=['accuracy'])
    h = model.fit(X_train, y_train_cat, batch_size=128*4, epochs=10, verbose=1, validation_split=0.3)

Train on 42000 samples, validate on 18000 samples
Epoch 1/10
42000/42000 [==============================] - 0s - loss: 0.5832 - acc: 0.8240 - val_loss: 0.4741 - val_acc: 0.8254
Epoch 2/10
42000/42000 [==============================] - 0s - loss: 0.2057 - acc: 0.9379 - val_loss: 0.1414 - val_acc: 0.9566
Epoch 3/10
42000/42000 [==============================] - 0s - loss: 0.1359 - acc: 0.9585 - val_loss: 0.1454 - val_acc: 0.9548
Epoch 4/10
42000/42000 [==============================] - 0s - loss: 0.0919 - acc: 0.9717 - val_loss: 0.0964 - val_acc: 0.9716
Epoch 5/10

gpu를 사용하면 아래와 같이 한 epoch 당 1초도 안되어서 트레이닝이되는 것을 볼 수 있습니다. 이번엔 cpu로 한 번 모델을 구축해보겠습니다. 아래와 같이 device 이름을 /cpu:0으로 변경한 후 모델을 구축, fitting을 하면, 한 epoch 당 2초가 걸리는 것을 확인할 수 있습니다.

import keras.backend.tensorflow_backend as K
with K.tf.device('/cpu:0'):
    model = Sequential()
    model.add(Dense(512, input_dim=28*28, activation='relu'))
    model.add(Dense(256, activation='relu'))
    model.add(Dense(128, activation='relu'))
    model.add(Dense(32, activation='relu'))
    model.add(Dense(10, activation='softmax'))

    model.compile(loss='categorical_crossentropy',
                  optimizer='rmsprop',
                  metrics=['accuracy'])
    h = model.fit(X_train, y_train_cat, batch_size=128*4, epochs=10, verbose=1, validation_split=0.3)

Train on 42000 samples, validate on 18000 samples
Epoch 1/10
42000/42000 [==============================] - 2s - loss: 0.6009 - acc: 0.8115 - val_loss: 0.3323 - val_acc: 0.8966
Epoch 2/10
42000/42000 [==============================] - 2s - loss: 0.2008 - acc: 0.9404 - val_loss: 0.3604 - val_acc: 0.8832
Epoch 3/10
42000/42000 [==============================] - 2s - loss: 0.1293 - acc: 0.9605 - val_loss: 0.1384 - val_acc: 0.9578
Epoch 4/10
42000/42000 [==============================] - 2s - loss: 0.0914 - acc: 0.9716 - val_loss: 0.1376 - val_acc: 0.9590
Epoch 5/10
42000/42000 [==============================] - 2s - loss: 0.0676 - acc: 0.9796 - val_loss: 0.1145 - val_acc: 0.9674
Epoch 6/10
42000/42000 [==============================] - 2s - loss: 0.0494 - acc: 0.9842 - val_loss: 0.0903 - val_acc: 0.9742

Python - 폴더 파일 리스트 가져오기

python에서 운영체제 폴더 아래의 파일들을 불러올 때는 os 패키지를 이용한다. 아래 코드에서 filepath 부분에 운영체제 상에 폴더 path를 적으면 폴더 내부의 파일 이름을 가져올 수 있다.

from os import listdir
from os.path import isfile, join
files = [f for f in listdir('filepath') if isfile(join('filepath', f))]

['Info.txt',
 'Licence.txt',
 'mdb001.pgm',
 'mdb002.pgm',
 'mdb003.pgm',
 'mdb004.pgm',
 'mdb005.pgm',
 'mdb006.pgm',
 'mdb007.pgm',
 'mdb008.pgm',

코드를 실행하면 폴더 내부의 파일들의 이름이 files 리스트에 저장된다. 만약 여기서 특정 문자열을 포함한 파일들만 남기고 싶으면 아래와 같이 문자열의 기본함수인 find함수를 이용한다. -1을 리턴하면 해당 문자열이 없는 것이기 때문에 아래의 if문은 해당 문자열을 포함한 원소만 가져오게 한다.

files = [x for x in files if x.find("mdb") != -1]
files

['mdb001.pgm',
 'mdb002.pgm',
 'mdb003.pgm',
 'mdb004.pgm',
 'mdb005.pgm',
 'mdb006.pgm',
 'mdb007.pgm',
 'mdb008.pgm'

Keras - CNN ImageDataGenerator 활용하기

keras에서는 이미지데이터 학습을 쉽게하도록 하기위해 다양한 패키지를 제공한다. 그 중 하나가 ImageDataGenerator 클래스이다. ImageDataGenerator 클래스를 통해 객체를 생성할 때 파라미터를 전달해주는 것을 통해 데이터의 전처리를 쉽게할 수 있고, 또 이 객체의 flow_from_directory 메소드를 활용하면 폴더 형태로된 데이터 구조를 바로 가져와서 사용할 수 있다. 이 과정은 매우 직관적이고 코드도 ImageDataGenerator를 사용하지 않는 방법에 비해 상당히 짧아진다. 환경은 keras tensorflow backend를 이용하였다.

1. 라이브러리 임포트

# Part 1 - Building the CNN

# Importing the Keras libraries and packages
from keras.models import Sequential
from keras.layers import Conv2D
from keras.layers import MaxPooling2D
from keras.layers import Flatten
from keras.layers import Dense

2. 모델 작성

ImageDataGenerator 클래스를 이용하는 것과 독립적으로 CNN 네트워크를 만들어준다. Input Image의 크기는 64x64이며, Conv-Pool-Conv-Pool-Fully Connected Layer로 이어지는 아주 간단한 구조의 CNN 모델이다. 데이터는 binary 데이터이기 때문에 activation function은 sigmoid함수 이며 손실 함수로 binary_crossentropy, 최적화 함수로는 adam을 사용한다.

# Initialising the CNN
classifier = Sequential()

# Step 1 - Convolution
classifier.add(Conv2D(32, (3, 3), input_shape = (64, 64, 3), activation = 'relu'))

# Step 2 - Pooling
classifier.add(MaxPooling2D(pool_size = (2, 2)))

# Adding a second convolutional layer
classifier.add(Conv2D(32, (3, 3), activation = 'relu'))
classifier.add(MaxPooling2D(pool_size = (2, 2)))

# Step 3 - Flattening
classifier.add(Flatten())

# Step 4 - Full connection
classifier.add(Dense(units = 128, activation = 'relu'))
classifier.add(Dense(units = 1, activation = 'sigmoid'))

# Compiling the CNN
classifier.compile(optimizer = 'adam', loss = 'binary_crossentropy', metrics = ['accuracy'])

3. ImageDataGenerator를 통한 트레이닝, 테스트 데이터 만들기

ImageDataGenerator를 만들 때 아래와 같이 rescale, shear_range, zoom_range, horizontal_flip을 설정해줌으로써 해당 데이터의 전처리를 어떻게 할지를 정해준다. 이 때 파라미터가 무엇을 뜻하는지를 정리하면 아래와 같다. (https://keras.io/preprocessing/image/)

• shear_range: Float. Shear Intensity (Shear angle in counter-clockwise direction as radians)
• zoom_range: Float or [lower, upper]. Range for random zoom. If a float,
• horizontal_flip: Boolean. Randomly flip inputs horizontally.
• rescale: rescaling factor. Defaults to None. If None or 0, no rescaling is applied, otherwise we multiply the data by the value provided (before applying any other transformation).

# Part 2 - Fitting the CNN to the images
from keras.preprocessing.image import ImageDataGenerator

train_datagen = ImageDataGenerator(rescale = 1./255,
                                   shear_range = 0.2,
                                   zoom_range = 0.2,
                                   horizontal_flip = True)

test_datagen = ImageDataGenerator(rescale = 1./255)

training_set = train_datagen.flow_from_directory('dataset/training_set',
                                                 target_size = (64, 64),
                                                 batch_size = 32,
                                                 class_mode = 'binary')

test_set = test_datagen.flow_from_directory('dataset/test_set',
                                            target_size = (64, 64),
                                            batch_size = 32,
                                            class_mode = 'binary')

classifier.fit_generator(training_set,
                         steps_per_epoch = 300,
                         epochs = 25,
                         validation_data = test_set,
                         validation_steps = 2000)

또한 flow_from_directory 메소드를 사용하면 폴더구조를 그대로 가져와서 ImageDataGenerator 객체의 실제 데이터를 채워준다. 이 데이터를 불러올 때 앞서 정의한 파라미터로 전처리를 한다. 마지막으로 fit_generator 함수를 실행함으로써 fitting이 이루어진다. 이 예제의 경우 ImageDataGenerator 객체를 사용하였기 때문에 fit_generator 함수를 통해 fitting한다. steps_per_epoch은 한 번 epoch 돌 때, 데이터를 몇 번 볼 것인가를 정해준다. [트레이닝데이터수/배치사이즈]를 사용하면 되는듯하다. validation_steps는 한 번 epoch 돌 고난 후, validation set을 통해 validation accuracy를 측정할 때 validation set을 몇 번 볼 것인지를 정해준다. 이 때도 마찬가지로 [validation data수/배치사이즈]를 사용하면 된다. 즉 ImageDataGenerator를 쓴 경우, fit_generator를 사용하면 된다고 기억하면 된다. 폴더 구조는 아래와 같이 하면 된다. flow_from_directory에 넣어준 경로(dataset/training_set) 밑에 이런식으로 class(cats, dogs) 별로 폴더를 만들고 폴더 밑에 이미지들을 넣어준다. 그러면 알아서 labeling을 하게 된다.

4. 테스트셋 예측하기

output = classifier.predict_generator(test_set, steps=5)
print(test_set.class_indices)
print(output)

{'cats': 0, 'dogs': 1}
[[ 0.58327454]
 [ 0.08756679]
 [ 0.44407782]
 [ 0.66460747]
 [ 0.38810217]
 [ 0.52190965]
 [ 0.67696238]
 [ 0.67470866]
 [ 0.77020335]
 [ 0.68121213]
 [ 0.66256768]
 [ 0.82018822]
...

Pandas 변수 정렬하기

Python에서 데이터 핸들링시 가장 많이 이용하는 Pandas 패키지를 이용하여 변수를 정렬하는 예제입니다. 예를 들어 아래와 같은 데이터셋이 있다고 합시다.

cluster  org      time
   1      a       8
   1      a       6
   2      h       34
   1      c       23
   2      d       74
   3      w       6 

import pandas as pd
data = pd.DataFrame({"cluster" : [1,1,2,1,2,3], "org":['a','a','h','c','d','w'], "time":[8,6,34,23, 74,6]})

import pandas as pd
data = pd.DataFrame({"cluster" : [1,1,2,1,2,3], "org":['a','a','h','c','d','w'], "time":[8,6,34,23, 74,6]})

이를 pandas DataFrame 객체로 읽기 위해서는 아래와 같은 구문으로 읽으면 됩니다.

import pandas as pd
data = pd.DataFrame({"cluster" : [1,1,2,1,2,3], "org":['a','a','h','c','d','w'], "time":[8,6,34,23, 74,6]})

이후에는 DataFrame 객체에 있는 sort_values를 호출하면 해당 변수에 대해 정렬을 할 수 있습니다. ascending 파라미터는 오름차순으로 정렬할지 여부를 결정합니다. 또한 여러개의 변수에 대해서 정렬하고 싶으면 list 자료구조에 변수를 넣어주기만 하면 됩니다.

data = data.sort_values(["time"], ascending=[False])

그러면 이와같이 결과가 나오는데 인덱스가 뒤죽박죽이 되어있다는 것을 알 수 있습니다. 기존 데이터셋의 인덱스를 그대로 갖고 있는 것입니다. 이로인해 plot을 할 때 문제가 생길 수 있으므로 인덱스를 초기화 하는 작업을 해줍니다. default는 drop=False인데 이렇게 하게 되면 기존의 인덱스가 하나의 컬럼으로 남게 됩니다. 많은 경우 기존의 인덱스는 필요없으므로 drop=True로 설정해줍니다.

data = data.reset_index(drop=True)