통계학

Quantile normalization

생물정보학에서 자주 사용되는 정규화 방법 중 하나인 quantile normalization이다. Quantile normalization은 비교하려는 샘플들의 분포를 완전히 동일하게 만들고 싶을 때, 또는 기준이 되는 분포(reference distribution)가 있는 경우 샘플들의 분포를 모두 기준 분포와 동일하게 만들고 싶을 때 사용할 수 있다.

복잡한 방식을 사용하는 과정이 아닌지라 쉽게 이해하고 쉽게 써먹을 수 있다.

 

알고리즘

  1. 샘플 내 값을 오름차순 정렬
  2. 아래 과정을 모든 n에 대하여 반복
    1. (기준 분포가 있을 떄)
      1. 모든 샘플의 n 번째 값들을 기준 분포의 n 번째 분위수로 대체
    2. (기준 분포가 없을 때)
      1. 각 샘플의 n 번째 값들을 평균(또는 중간값)
      2. 각 샘플의 n 번째 값들을 모두 2-1.의 값으로 대체
  3. 정렬된 데이터를 원래 데이터의 순서대로 복구

 

이해가 어렵다면 바로 아래 “예시”로 넘어가도 된다.

“값들을 순서대로 나열했을 때 n 번째 값”이 바로 통계학에서의. n 번째 분위수(n-th quantile)의 개념에 해당하고, 각 샘플 분포의 n 번째 분위수를 동일하게 조정해주는 작업이므로 quantile normalization이라는 이름이 붙었다.

딱히 알고리즘이랄 것도 없다. 1.에서 정렬 시 오름차순 대신 내림차순을 써도 아무 상관이 없다. 기준 분포가 없을 때 2-1.에서 평균, 중간값 대신 가중평균 등 자신에게 적합한 대푯값을 사용해도 된다.

 

예시

1. 이해를 위한 간단한 예시

다음 테이블과 같은 세 샘플(A, B, C)의 데이터가 있다고 하자. A, B, C의 분포가 동일해지도록 median quantile normalization해보자. 셀 색상은 원래 데이터의 순서를 표시하기 위해 넣어보았다.

qn_whole

  1. 열 별로 값을 정렬한다.
  2. 정렬된 데이터에서 행 별로 중간값을 계산한다.
  3. 각 행의 값을 2.에서 계산한 중간값으로 바꾼다.
  4. 정렬을 해제하여 기존 데이터와 같은 순서로 돌려준다

 

2. 규모가 좀 더 큰 예시

서로 다른 샘플에 어떤 실험을 한 결과 데이터가 세 세트(A, B, C) 있다고 하자. 세 샘플을 용이하게 비교하기 위해 샘플의 분포를 quantile normalize하여 동일하게 만들어보자.

 

기준 분포가 있을 때

표준정규분포를 기준 분포로 삼는다고 하자.

정규화 전 각 샘플의 분포는 아래와 같다.

A  B   C
0   3.570614    5.015660    6.526778
1   5.965098    1.258324    2.285988
2   3.767857    0.974910    4.374467
3   3.538398    3.362349    1.779624
4   5.178286    1.575602    3.065568
5   3.578033    3.434900    5.086503
6   5.330161    3.431059    3.573926
7   5.002079    4.885904    2.318780
8   5.694074    4.000688    4.615286
9   6.026262    -0.499165   1.999979
before

 

표준정규분포의 101-quantile을 사용하여 quantile normalization하면,

A  B   C
0   -1.000990   1.410420    1.410420
1   0.960838    -1.180947   -1.000990
2   -0.922178   -1.346263   0.162024
3   -1.086568   0.162024    -1.232341
4   0.162024    -0.848716   -0.476577
5   -0.960838   0.263612    0.620918
6   0.289397    0.238000    -0.187224
7   -0.012409   1.232341    -0.922178
8   0.651302    0.682300    0.368003
9   1.042824    -2.330079   -1.132497
after_ref

모든 샘플의 분포가 표준정규분포와 일치하는 것을 확인할 수 있다.

 

기준 분포가 없을 때

기준 분포가 없이 median을 사용하여 정규화하면,

A  B   C
0   -1.000990   1.410420    1.410420
1   0.960838    -1.180947   -1.000990
2   -0.922178   -1.346263   0.162024
3   -1.086568   0.162024    -1.232341
4   0.162024    -0.848716   -0.476577
5   -0.960838   0.263612    0.620918
6   0.289397    0.238000    -0.187224
7   -0.012409   1.232341    -0.922178
8   0.651302    0.682300    0.368003
9   1.042824    -2.330079   -1.132497
after

역시 마찬가지로 모든 샘플의 분포를 동일하게 만들 수 있다.

 

참고

이 코드를 돌리면 예시 2-2를 직접 해볼 수 있다.

Sequence Motif 등장 확률

전체 n ntd로 이루어진 target RNA 한 가닥이 있다고 하자. 이 RNA는 완전히 랜덤하게 만들어진 가닥이라고 가정한다. 즉, 각 위치에 A, U, G, C가 같은 확률로 존재할 수 있다.

이제 k-mer sequence motif가 하나 있다고 하자. 이 sequence motif는 특정되어 있다.

이 때 Target RNA 가닥에서 이 motif가 한 번이라도 등장할 확률은 얼마일까?

우연한 기회에 이 문제에 대해 생각해보게 되어 여사건의 확률 등의 방법을 생각해보다 이런 방법으로는 쉽게 풀리지 않는다는걸 알게 되었다.

검색을 해봐도 정확한 해답을 제시하는 사람은 많지 않았다 (너무 쉬운 문제여서 인지?). 여기서는 확률의 inclusion-exclusion principle을 이용한 해답을 소개하려 한다.

 

마주서기 문제

본 문제와 비슷하지만 비교적 단순한 문제를 먼저 살펴보자.

n쌍의 부부가 있다. 한 쪽에는 남편끼리, 다른 쪽에는 부인끼리 늘어 설 때, 적어도 한 쌍의 부부가 서로 마주보고 설 확률은 얼마일까?

이 문제 역시 여사건의 확률을 이용해서 풀기는 쉽지 않다 (재귀함수를 만들게 된다). 대신 inclusion-exclusion principle을 사용해서 풀어보자.

사건 A_i를 “i번째 부부가 서로 마주보고 서는 사건”이라고 하자. 그러면

\begin{array}{lcl} &&P(A_i) = \dfrac{(n-1)!}{n!} \text{, } i=1, ..., n \\ &&P(A_i \cap A_j) = \dfrac{(n-2)!}{n!} \text{, } i < j \\ && ... \\ &&P(A_1 \cap A_2 \cap ... \cap A_n) = \dfrac{1}{n!} \end{array}

이므로

\begin{array}{lcl} P(\cup_{i=1}^{n} A_i) &=& \sum_{i=1}^{n}P(A_i) - \sum_{i<j}P(A_i \cap A_j) \\ & & + ... + (-1)^{n+1}\sum_{i<j<...<l}P(A_i \cap A_j \cap ... \cap A_l) \\ &=& \sum_{i=1}^{n} (-1)^{i+1} \dfrac{1}{i!} \end{array}

이다.

 

본 문제

사건 A_i \text{, } i=1, ..., m를 “i번째 위치에 motif가 등장하는 사건”이라고 하자. 단 mnk로 나눈 몫이다. 그러면

\begin{array}{lcl} &&P(A_i) = \dfrac{4^{n-k}}{4^n} \text{, } i \leq n-(k-1) \\ &&P(A_i \cap A_j) = \dfrac{4^{n-2k}}{4^n} \text{, } i<j \text{ and } i, j \leq n-2(k-1) \\ &&P(A_i \cap A_j \cap A_k) = \dfrac{4^{n-3k}}{4^n} \text{, } i<j<k\text{ and } i, j, k \leq n-3(k-1) \\ &&... \\ &&P(A_1 \cap A_2 \cap ... \cap A_m) = \dfrac{4^{n-mk}}{4^n} \end{array}

이므로 마주서기 문제에서와 같은 방법으로

P(\cup_{i=1}^{m} A_i) = \sum_{i=1}^{m} (-1)^{i+1} \binom{n-i(k-1)}{i} 4^{-ki}

임을 알 수 있다.

 

이 방법으로 1~500 ntd의 random target sequence에 3, 4, 7-mer fixed motif가 존재할 확률을 그래프로 나타내보면 다음과 같다.

Thu Oct 25 00:15:46 2018

 


참고

 

 

Batch Normalization 이해하기

현대적인 딥러닝 모델을 디자인할 거의 항상 빠지지 않고 쓰이는 테크닉들이 있다. 하나는 recurrent 구조 (LSTM, Attention)이고 다른 하나는 batch normalization (BatchNorm)이다. LSTM과 attention 대해서는 recurrent neural net 다루면서 자세히 살펴보도록 하고 이번 글에서는 학습 과정에서 뉴럴넷을 안정시켜주는 표준화 기법 하나인 batch normalization 대해 다뤄보겠다.

 

  • 기존 방법의 문제점
  • BatchNorm
    • 알고리즘
    • 테스트할 때
    • BN layer
  • TensorFlow 구현

 

기존 방법의 문제점

BatchNorm이 어떤 의미를 가지는지를 알기 위해서는 BatchNorm이 고안되기 이전의 딥러닝 모형 초기화 및 학습 과정 표준화 과정을 둘러볼 필요가 있다.

뉴럴넷이 안정적으로 잘 학습되기 위해서는 입력층에 넣을 인풋과 각 층의 weight를 표준화할 필요가 있다. BatchNorm이 고안되기 전에는 두 가지 방법을 주로 사용했는데, 이전 포스트[1, 2]에서 각각의 방법을 간단히 다룬 바 있다. 간단히 복기하자면 이렇다: (1) 인풋은 centering scaling하고 (2) 인풋 뉴런 n개인 층의 weight \div \sqrt{n/2}로 표준화한다. 단순한 방법이지만 표준화하지 않은 입력, 가중치값을 사용했을 때에 비해 더 빨리, 더 좋은 성능으로 수렴하는 것을 경험적으로 확인할 수 있다.

여기서 중요한 문제가 발생한다. 입력층에 넣는 인풋은 표준화할 수 있다. 뉴럴넷에 넣기 전에 우리가 원하는 방식으로 원하는 만큼 preprocessing을 하면 된다. 그 결과 입력층의 input distribution은 항상 비슷한 형태로 유지가 되고 안정적으로 가중치 학습을 진행할 수 있다.

e18489e185b3e1848fe185a6e1848ee185b53.png

그러나 은닉층은 인풋의 분포가 학습이 진행됨에 따라 계속 변한다. 은닉층은 이전 레이어의 activation f(XW)을 입력으로 받는다. 학습 과정에서 가중치 W의 값이 W^\prime로 업데이트되면 이전 레이어의 activation 또한 f(XW^\prime)로 바뀌게 된다. 은닉층의 입장에서는 인풋 값의 분포가 계속 널뛰는 것이나 마찬가지이다. 입력 분포의 형태가 유지되지 않으므로 학습도 잘 진행되지 않는다. 그라디언트 값이 큰 학습 초기일수록 문제가 더 심각해진다.

스케치

 

Batch Normalization

알고리즘

바로 위에서 언급한 문제를 internal covariate shift라고 한다. 그대로 입력층보다 깊은, 내부에 있는(internal) 층의 입력값, 공변량(covariate) 고정된 분포를 갖지 않고 이리저리 움직인다(shift) 의미이다. BatchNorm 바로 internal covariate shift 해결하는 테크닉이다.

[1]

은닉층의 입력도 표준화한다면 안정적으로 깊은 레이어의 가중치도 학습시킬 수 있을 것이다. “은닉층의 입력을 표준화한다는 것은 곧이전 층의 출력(raw activation)을 표준화한다는 의미와 같다.

딥러닝은 거의 항상 전체 샘플을 mini batch로 나누어 학습하고 가중치를 업데이트하므로 이전 층의 raw activation을 표준화할때도 각 batch마다 따로 표준화하면 된다.

스케치

이와 같이 각각의 minibatch 평균 \mu_{\mathcal{B}} = \frac{1}{m} \sum_i {x_iw_i} 표준편차 \sigma_{\mathcal{B}} = \frac{1}{m} \sum_i {(x_iw_i - \mu_{\mathcal{B}})^2} 표준화한 activation a_s = f(\frac{XW_1 - \mu_{\mathcal{B}}}{\sigma_{\mathcal{B}}}) 은닉층 B 입력으로 사용하면 은닉층 B 입력은 고정된 분포를 따른다.

쉬워도 너무 쉽다. 이렇게만 하면 될 것 같지만..

[1 문제점]

문제가 가지 있다. 이렇게 은닉층의 입력을 표준화하면 gradient update 과정에서 bias(편향)값이 무시된다. [1]만을 사용해서 표준화한다고 그라디언트 업데이트 과정을 자세히 살펴보자. Raw activation a_r = wx + b라고 E(a_r) = \frac{1}{n} \sum_i a_{r_{i}}이므로

  1. 그라디언트를 계산한다.
    • \Delta b \propto - {\partial L}/{\partial b},  where L is a loss function.
  2. 편향(과 가중치)을 업데이트한다.
    • b \gets b + \Delta b
  3. 편향을 업데이트한 이후의 raw activation:
    • a_r ^\prime = wx + (b + \Delta b)
  4. [1] 이용해서 센터링만 raw activation:
    • \begin{array}{lcl} a_{r_{centered}} ^\prime &=& a_r ^\prime - E(a_r ^\prime) \\ &=& \{(wx + b) + \Delta b\} - \{ E[wx + b] + \Delta b \} \\ &=& (wx + b) - E[wx + b] \end{array}

Bias b 업데이트 \Delta b 완벽하게 캔슬되었다. 초기 편향값에서 이상 업데이트가 되지 않는 것이다. 종류의 파라미터 w, b 사용했는데 파라미터 w 가지만 사용하는 단순한 모형으로 irreversible하게 변환된 것이다.

이 때문에 b 대신 편향의 역할을 할 파라미터를 추가해야한다. 이 파라미터는 그라디언트 업데이트 과정에서 무시되어서는 안된다.

다른 문제도 있다. raw activation 분포를 고정시키는 것은 좋지만 항상 N(0, 1) 고정시킬 필요는 없다. 적절하게 scaling, shifting activation \gamma \cdot \frac{a_r - \mu_{\mathcal{B}}}{\sigma_{\mathcal{B}}} + \beta 사용하는 것이 학습에 도움될 수도 있다.

형태의 activation 사용할 경우 필요하다면 표준화를 되돌릴 수도 있다. \gamma = \sigma_{\mathcal{B}}, \beta = \mu_{\mathcal{B}} \gamma \cdot \frac{a_r - \mu_{\mathcal{B}}}{\sigma_{\mathcal{B}}} + \beta = a_r이기 때문이다.

[2]

위의 문제를 극복하기 위해 표준화한 scaling shifting raw activation, 즉

a_{BN} = \gamma \cdot \frac{XW_1 - \mu_{\mathcal{B}}}{\sigma_{\mathcal{B}}} + \beta

activation function f 입력으로 사용한다. 은닉층 B 입력으로는 f(a_{BN}) 사용한다. 방법을 BatchNorm이라고 한다. \gamma, \beta 파라미터로 학습 과정에서 업데이트되는 값이다.

BatchNorm 장점이 꽤나 많은데

  • bias 업데이트를 무시하지 않는다. \beta bias처럼 행동한다. \beta 업데이트는 표준화해도 캔슬되지 않는다.
  • 은닉층마다 적절한 input distribution 가질 있다. scaling factor \gamma shifting factor \beta 사용해서 적절한 모양으로 입력분포를 조정할 있다.
  • 필요한 경우 표준화를 하지 않을 수도 있다. 위에서 언급한 \gamma = \sigma_{\mathcal{B}}, \beta = \mu_{\mathcal{B}} 경우이다.
  • Activation 값을 적당한 크기로 유지하기 때문에 vanishing gradient 현상을 어느정도 막아준다. 덕분에 tanh, softmax같은 saturating nonlinearity 사용해도 문제가 생긴다.
  • batch-wise로 계산하기 때문에 컴퓨팅하기 용이하다.
  • 위의 장점들을 모두 가지면서, 동시에 층마다 입력 분포를 특정 형태로 안정시켜서 internal covariate shift 방지할 있다.
  • 입력 분포가 안정되므로 학습시 손실함수가 더 빨리, 더 좋은 값으로 수렴한다.
  • 초기 learning rate를 크게 설정해도 안정적으로 수렴한다고 한다.
  • Weak regularizer로도 작용한다고 한다.

이쯤 되면 거의 만능이다.

테스트

지금까지 다룬 내용은 모두 학습 과정에서 일어나는 일들이다. 학습 과정에서는 raw activation minibatch mean, stdev 표준화하면 됐었다. 그런데 학습을 마치고 테스트(또는 evaluation, inference) 때에는 minibatch mean, stdev 존재하지 않는다.

테스트 과정에서는 대신 전체 training data mean, stdev 사용해서 BatchNorm 한다. 전체 training data mean, stdev 번에 계산하기에는 메모리의 제약이 있으므로, minibatch statistic 평균낸 값을 대신 사용한다.

, n개의 minibatch 있을 ,

\hat{\mu} = \frac{1}{n} \sum_i {\mu_{\mathcal{B}}^{(i)}}
\hat{\sigma} = \frac{1}{n} \sum_i {\sigma_{\mathcal{B}}^{(i)}}

Minibatch statistic 따로 저장할 필요 없이 학습 과정에서 moving average \hat{\mu}, \hat{\sigma} 계산하면 된다. Exponential moving average 사용해도 좋다.

i번째 minibatch statistic 각각 \mu_{\mathcal{B}}^{(i)}, \sigma_{\mathcal{B}}^{(i)}라고 ,

\hat{\mu} \gets \alpha \hat{\mu} + (1-\alpha) \mu_{\mathcal{B}}^{(i)}
\hat{\sigma} \gets \alpha \hat{\sigma} + (1-\alpha) \sigma_{\mathcal{B}}^{(i)}

BatchNorm layer

ReLU activation 뉴럴넷의 레이어로 나타낼 있듯 BatchNorm 또한 레이어로 표현할 있다. BN layer raw activation activation function 사이에 위치한다. Convolutional layer에 BatchNorm을 적용하고 싶을 때에도 동일하게 raw feature map과 ReLU layer 사이에 BN layer를 추가하면 된다.

e18489e185b3e1848fe185a6e1848ee185b57.png

BN layer mini batch raw activations a_r 입력받아 아래와 같은 연산을 수행하여 다음 레이어(activation function f) 전달한다.

BN_{\gamma, \beta}(a_r) = \gamma \cdot \frac{a_r - \mu_{\mathcal{B}}}{\sigma_{\mathcal{B}}} + \beta

또한 테스트 사용하기 위해 학습 과정에서 minibatch statistic exponential moving average(또는 그냥 MA) minibatch마다 업데이트한다.

 

TensorFlow 구현

구글에서 고안한 방법답게 TensorFlow에 이 내용들이 친절히 함수로 구현되어 있다. tf.nn.batch_normalization, tf.contrib.slim.batch_norm를 쓰면 간단히 위 알고리즘을 모형 구축에 사용할 수 있다.

tf.nn.batch_normalization을 사용할 경우, minibatch statistic의 EMA를 계산하는 코드를 따로 작성해야 한다.

tf.contrib.slim.batch_norm를 사용할 경우 is_training 옵션을 True로 주면 자동으로 EMA를 계산해서 저장하고, False로 주면 저장된 EMA 값으로 activation을 표준화한다.

TF-Slim 레이어에도 쉽게 적용시킬 수 있다.

import tensorflow as tf
import tensorflow.contrib.slim as slim

bn_params = {"decay": .9,
             "updates_collections": None,
             "is_training": tf.placeholder(tf.bool)}
net = slim.fully_connected(input, 1024,
                           normalizer_fn=slim.batch_norm,
                           normalizer_params=bn_params)

Convolutional layer에도 마찬가지다.

net = slim.conv2d(input, 64, [5,5], padding="SAME",
                  normalizer_fn=slim.batch_norm,
                  normalizer_params=bn_params)

 

 

참고

Logistic, cross-entropy loss의 확률론적 의미

머신러닝을 공부하다보면 logistic을 참 많이 접하게 된다. 당장 로지스틱 회귀만 해도 그렇고, 딥러닝에서 자주 사용하는 saturating non-linearity 중 하나인 softmax function도 logistic function의 multi-class 버전이니 말이다.

Logistic(또는 softmax) function은 어떤 값들을 (0, 1) 사이의 값으로 squash하면서, 동시에 값들의 총 합이 1이 되도록 하는 성질이 있는데, 이 결과 logistic function의 함수값이 Kolmogorov의 확률의 정의에 부합하기 때문에 주로 값들을 확률화하는 함수로 쓰인다. 딥러닝에서 분류문제를 풀때 output layer의 activation function으로 자주 쓰이곤 하는 이유도 여기에 있다 (각 클래스의 ‘확률’ 값을 계산하기 위해).

그런데 하고많은, (0, 1)에 bounded되어있고 합이 1이 되는 함수들 중에서 굳이 logistic function을 사용하는 이유는 뭘까? 그건 바로 logistic function이 확률값에서부터 자연스럽게 유도되는, 내재적으로 확률의 성질을 가지는 함수이기 때문이다. 이 글은 logistic이 그냥 이유없이 확률처럼 생겨서 만들어지고 사용되는 함수가 아니라는 것을 증명해보는 지루한 글이다. Logistic function은 내재적으로 cross-entropy loss와 깊은 연관이 있는데, 이 로스를 줄이는 것이 왜 classification을 해결하는 과정이 되는지 역시 같은 맥락에서 함께 살펴볼 것이다.

 

Odds, Log odds, Logit, 그리고 Logistic

어떤 사건이 일어날 확률은 P(x) 외에도 다양한 방법으로 표현될 수 있다. Odds가 그 종류 중 하나다. 주로 질병통계학이나 도박과 같이 카테고리가 두 개(성공과 실패/발병과 비발병)인 분야에서 자주 사용하는 표현인데, 로지스틱 함수의 확률론적 의미는 바로 여기에서부터 시작한다.

발생 확률(또는 성공 확률)이 p인 어떤 사건 A의 Odds는 다음과 같이 나타낼 수 있다.

Odds(A) = \frac{p}{1-p}

즉, Odds는 성공 확률이 실패 확률에 비해 얼마나(몇 배나) 더 큰가?를 나타내는 척도라고 해석할 수 있다. Odds가 클수록 성공확률이 큰 것이다. Odds는 평범한(?) 확률값 p에 비해 성공확률값의 변화에 민감하다는 특징이 있다. 또한 질병통계학에서는 relative risk 등의 다른 표현에 비해 robust하기에 자주 사용한다.

Odds에 (자연)로그를 취한 것을 Log odds라고 한다.

\log{Odds(A)} = \log{\frac{p}{1-p}} = \log{p} - \log{(1-p)}

로그변환은 통계학에서 자주 사용하는 변환으로,

  • 함수의 증감 형태, convex/concave 형태를 유지시키고
  • 극점의 위치를 유지시키며
  • 곱(또는 나눗셈)으로 표현된 식을 선형조합의 꼴로 풀어쓸 수 있도록 해준다

는 장점이 있다. 즉, 계산은 용이하게 해주는데 함수의 특성은 그대로 유지하는 변환이라 할 수 있겠다. 로그변환의 이러한 성질에 의해 Log odds는 여전히 값이 클수록 성공확률이 큰 것을 의미하는 함수이다.

한 편, Logit function은 다음과 같은 꼴의 함수를 말한다. Logit function은 (0, 1)을 support로 가지며 (- \infty, + \infty) 값을 domain으로 가진다.

logit(x) = \log{\frac{x}{1-x}}

어디서 많이 본 꼴이다. 그렇다. Log odds는 성공확률 p에 대한 Logit function 함수값과 일치한다 (\log{Odds(A)} = logit(p)). 다시말해 Logit 함수는 log Odds의 성질을 가지는 함수이다.

Logistic function은 바로 이 Logit function의 역함수이다.

 

logistic(x) = \frac{1}{1+e^{-x}}

성공확률 p를 갖는 사건 A를 바탕으로 해석해보면 Logistic function은 성공확률 p 그 자체이다.

\frac{p}{1-p} = e^{\log{Odds(A)}},

p = e^{\log{Odds(A)}} -e^{\log{Odds(A)}}p,

(1+e^{\log{Odds(A)}})p = e^{\log{Odds(A)}},

\begin{array}{lcl}\therefore p &=& \frac{e^{\log{Odds(A)}}}{1+e^{\log{Odds(A)}}} \\ &=& \frac{1}{1+e^{-\log{Odds(A)}}} \end{array}

 

Binary classification과 Logistic function

두 개의 클래스(클래스 1, 2)를 분류하는 문제를 생각해보자. 우리가 분류에 사용할 데이터를 X,특정 샘플이 클래스 1, 2에 속한다고 판단하는 사건을 각각 Y_1, Y_2라고 하자.

Binary classification 문제는 데이터 X가 주어졌을 때 둘 중 어느 클래스에 속할지 그 확률 P(Y_i|X) (posterior probability)를 구하는 문제라고 할 수 있다. P(Y_1|X) >P(Y_2|X)일 때 우리는 Y_1이 옳은 클래스일 것이라고 예측을 하게 되고, P(Y_1|X) < P(Y_2|X)일 때 Y_2가 옳은 클래스일 것이라고 예측하게 된다.

Bayes’ rule에 의해서 posterior probability P(Y_i|X)는 likelihood P(X|Y_i)와 prior probability P(Y_i)의 곱에 비례한다. (분모인 prior predictive probability P(X)는 클래스에 관계 없이 X의 marginal probablity로 항상 같으므로 덜 중요하다)

P(Y_i|X) = \begin{cases}P(Y_1|X)=\frac{P(X|Y_1)P(Y_1)}{P(X)} =\frac{P(X|Y_1)P(Y_1)}{P(X|Y_1)P(Y_1) +P(X|Y_2)P(Y_2)}  \\ P(Y_2|X)=\frac{P(X|Y_2)P(Y_2)}{P(X)} =\frac{P(X|Y_2)P(Y_2)}{P(X|Y_1)P(Y_1) +P(X|Y_2)P(Y_2)} \end{cases}

따라서

P(Y_i|X) = \begin{cases}P(Y_1|X) \propto P(X|Y_1)P(Y_1)  \\ P(Y_2|X) \propto P(X|Y_2)P(Y_2) \end{cases}

여기서 a_i = \log{P(X|Y_i)P(Y_i)}를 정의하면 다음과 같이 각 클래스의 posterior probability를 표현할 수 있게 된다.

\begin{array}{lcl} P(Y_1|X) &=& \frac{e^{a_1}}{e^{a_1} + e^{a_2}} \\ &=& \frac{1}{1+e^{-(a_1-a_2)}} \end{array}

a = a_1 - a_2 = \log{ \frac{P(X|Y_1)P(Y_1)}{P(X|Y_2)P(Y_2)} }라고 하면 P(Y_1|X)a에 대한 로지스틱 함수값이 된다.

P(Y_1|X) =\frac{1}{1+e^{-a}} = logistic(a)

그런데 a는 위의 정의에 따라

\begin{array}{lcl} a &=& \log{ \frac{P(X|Y_1)P(Y_1)}{P(X|Y_2)P(Y_2)} } \\ &=& \log{ \frac{P(Y_1|X)}{P(Y_2|X)} } \\ &=& \log{ \frac{P(Y_1|X)}{1-P(Y_1|X)} } \\ &=& \log{Odds(Y_1|X)} \end{array}

즉, 데이터 X를 바탕으로 Y_1이 옳은 클래스라고 판단할 log odds이다.

위에서 사건 A의 logistic은 사건의 성공확률 p 그 자체라는 것을 보였으므로, 결론적으로 logistic(Y_1|X)는 사건 Y_1|X의 확률, 즉 내재적으로 P(Y_1|X) 그 자체이다 (바로 위의 증명과 완전히 일맥상통한다).

이때문에 우리는 확률분포에 대한 가정 없이도 (deterministic하게) logistic function을 각 클래스로 판단할 ‘확률값’으로 곧바로 사용할 수 있는 것이다.

Multi-class classification과 Softmax function

n 개의 클래스(클래스 1, 2, … n)를 분류하는 문제를 생각해보자. 우리가 분류에 사용할 데이터를 X,특정 샘플이 클래스 1, 2, … n에 속한다고 판단하는 사건을 각각 Y_1, Y_2, ..., Y_n라고 하자. Logistic function을 유도했을 때와 마찬가지로,

P(Y_i|X) = \begin{cases}P(Y_1|X)=\frac{P(X|Y_1)P(Y_1)}{P(X)} =\frac{P(X|Y_1)P(Y_1)}{\sum_{i=1}^{n} P(X|Y_i)P(Y_i)}  \\ P(Y_2|X)=\frac{P(X|Y_2)P(Y_2)}{P(X)} =\frac{P(X|Y_2)P(Y_2)}{\sum_{i=1}^{n} P(X|Y_i)P(Y_i)} \\ ... \\ P(Y_n|X)=\frac{P(X|Y_n)P(Y_n)}{P(X)} =\frac{P(X|Y_n)P(Y_n)}{\sum_{i=1}^{n} P(X|Y_i)P(Y_i)} \end{cases}

이고, a_i = \log{P(X|Y_i)P(Y_i)}를 정의하면 posterior probability는 다음과 같이 Softmax function의 꼴로 표현된다.

\begin{array}{lcl} P(Y_1|X) &=& \frac{e^{a_1}}{e^{a_1} + e^{a_2} + ... +e^{a_n}} \\ &=&\frac{e^{a_1}}{\sum_{i=1}^{n} e^{a_i}} \end{array}

이처럼 Logistic function을 multi-class 문제로 일반화시키면 Softmax function을 얻을 수 있다. 이때문에 Softmax function을 multi-class logistic function이라고 하기도 한다.

 

Cross-entropy loss와 MLE

i번째 관측치 벡터 (x_i)의 ground truth 클래스를 t_i, 분류기를 통해 판단한 클래스를 y_i라고 하면,

t_i|x_i, a {\sim}^{iid} Bernoulli(p_i),  where p_i = P(t_i=1|x_i) = P(Y_1|X) = y_i

라고 가정할 수 있다. 이 때 likelihood는 다음과 같이 적을 수 있다.

\begin{array}{lcl} L = \prod_{i=1}^{n} P(t_i|x_i, a) &=& \prod_{i=1}^{n} p_{i}^{t_i} (1-p_i)^{1-t_i} \\ &=& \prod_{i=1}^{n} y_{i}^{t_i} (1-y_i)^{1-t_i}  \end{array}

Maximum likelihood estimation(MLE)을 위해서 log likelihood \log{L}을 구하면

\log{L} = \sum_{i=1}^{n} \log{\{y_{i}^{t_i} (1-y_i)^{1-t_i}\}} = \sum_{i=1}^{n} \{ t_i \log{y_i} + (1-t_i) \log{(1-y_i)} \}

인데, 이 식에 negative(-)를 취하면 cross-entropy loss(CE)의 정의와 일치한다. 따라서

\hat{\theta}^{MLE} = argmax_{\theta}(\log{L}) = argmin_{\theta}(-\log{L}) =argmin_{\theta}(CE)

이므로, 위의 Bernoulli 분포 가정 하에서 cross-entropy loss를 minimize하는 것은 MLE와 동일한 과정이 된다. 이때문에 cross-entropy loss를 최소화하는 과정을 통해 올바른 클래스 t_i에 가까운 예측 확률값 y_i를 얻을 수 있게 된다.

 

참고

  • 이 글을 많이 참고했다. 이 글을 바탕으로 노테이션을 정리하고 수식 유도를 풀어 쓴 정도다.
  • Odds에 대해서는 이 글을 참고했다.
  • 이 글에서 영감(?)을 얻었다.

EFA로 내게 맞는 스마트폰 카메라 찾기

해를 거듭하면서 스마트폰 성능이 눈에 띄게 향상되고 있다. 특히나 아이폰에 탑재되는 A-시리즈 APU의 성능은 작년부터 PC를 위협해오더니 급기야 13인치 맥북 프로를 벤치마크 상에서 앞지르기에 이르렀다. 일정 가격 이상의 대부분의 스마트폰 성능이 상향 평준화 되어버린 지금, Kantar가 2014년 수행한 북미 시장에서의 설문조사를 보면 스마트폰 구매 선택은 성능 외적인 부분에서 주로 결정되는 것으로 보인다.

특히 카메라 성능은 고가의 스마트폰을 구입하는 데에 있어서 결정적인 역할을 하는 요인으로 생각된다. 2016년 이전에는 코빼기도 보이지 않던 DxOMark 카메라 벤치마크가 이제는 많은 리뷰에서 빠짐없이 등장하는 현상이 이를 뒷받침한다. 비록 DxOMark에서 다양한 성능 카테고리에 대해서 카메라를 테스트하기는 하지만, 벤치마크와 실제 사용자 경험은 차이가 있을 수 있다. 테스트 점수가 높다고해서 내게 만족스러운 카메라가 되는 것은 아니라는 의미이다.

어떤 카메라가 내게 맞는 카메라인지에 대해 대략의 답을 제공하기 위해 DxOMark 벤치마크 데이터에 대해 탐색적 요인 분석(exploratory factor analysis; EFA)을 해보았다.

분석을 시작하면서 답을 찾고자 한 질문들은 이렇다.

  1. 각 브랜드는 스마트폰 카메라를 제작할 때 어떤 성능요인(색감, 선명도, …)에 중점을 두는가
  2. 각 브랜드의 스마트폰 카메라는 시간에 따라 어떤 방향으로 발전해왔는가
  3. 어떤 스마트폰 브랜드끼리 서로 비슷한 사진을 찍어내는가

주로 애플과 삼성을 중점으로 분석을 진행했다. 이 외의 브랜드의 스마트폰 카메라 분석을 원한다면 다음 gist를 참조하시길. 분석과 코드를 정리해 두었다 [벤치마크 크롤러][노트북].

 

최근 데이터들의 3-factor EFA 결과

DxOMark는 2017년 9월에 새로운 테스트 프로토콜을 제시한 바 있다. 테스트 프로토콜에 따른 테스트 점수 변동이 있을 수 있으므로, 이전 프로토콜에서 만들어진 데이터(이 글에서는 “과거” 데이터라 호칭한다)와 새로운 프로토콜(“최근”)에서 만들어진 데이터를 따로 분석하기로 했다.

Varimax 방법으로 새 프로토콜에서 얻은 데이터들을 회전시킨 EFA 결과, 다음과 같은 factor loading을 얻을 수 있었다.

output_39_0

각 요인에 크게 연관되어 있는 manifest variable에 따라 요인을 해석해서 이름을 붙였다.
FA1: AF, texture, exposure&contrast 등에 강한 음의 상관관계가 있으므로 “Bad Focus”.
FA2: artifacts, noise 등에 강한 양의 상관관계가 있으므로 “Sharp, Clear”.
    FA3: flash, texture 등에 강한 양의 상관관계가 있으므로 “High Detail”.
즉, FA1 값이 클수록 포커스가 잘 안잡히는 카메라, FA2가 클 수록 선명하고 또렷한 카메라다.

현재 시점 (2017/12/29)에 DxOMark에 올라온 모든 최근 스마트폰 카메라의 점수를 요인 점수로 변환해서 그래프로 나타내 보았다.

output_43_0

같은 색 포인트는 같은 회사에서 제작된 스마트폰이다. 이 그래프에서 각 브랜드의 장단점을 파악할 수 있다. 예를들어 구글의 경우 Sharp, Clear 점수가 조금 떨어지지만 포커싱과 디테일이 매우 우수하다. 반면 애플은 포커싱이 비교적 안좋은 반면 Sharp, Clear과 디테일이 우수한 것을 볼 수 있다. 삼성과 구글의 지향점이 서로 비슷한 듯하다.

각 브랜드의 최근 카메라 제작 트렌드를 좀 더 자세히 들여다보기 위해 삼성과 애플 제품만 따로 나타내보았다.

output_42_0

애플 아이폰 모델이 최신형에 이르면서 포커싱보다는 Sharp, Clear과 디테일에 중점을 두고 발전하는 것을 확인할 수 있다. 삼성 또한 갤럭시 S6 엣지에서 갤럭시 노트 8에 이르면서 Sharp, Clear보다는 포커싱과 디테일에 중점을 두는 모양새를 볼 수 있다.

 

과거 데이터들의 3-factor EFA 결과

최근 데이터와 동일한 방법으로 분석을 진행했다. Varimax rotated factor loading은 아래와 같다.

output_51_0

최신 데이터와는 조금 다른 해석을 사용해서 요인에 이름을 붙였다.
FA1: AF, texture, noise 등에 강한 음의 상관관계가 있으므로 “Bad Focus”.
FA2: color, exposure&contrast 등에 강한 양의 상관관계가 있으므로 “Colorful”.
    FA3: noise, flash, exposure&contrast 등에 강한 양의 상관관계가 있으므로 “Clear”.
즉, FA1 값이 클수록 포커스가 잘 안잡히는 카메라, FA2가 클 수록 선명하고 또렷한 카메라다.

DxOMark에서 테스트한 과거 모든 스마트폰 카메라를 스캐터플롯으로 나타내어 보았다.

output_54_0

삼성과 애플 제품들만 모아보았다.

output_55_0

삼성 제품들과 애플 제품들이 서로 다른 클러스터로 구별된다.
삼성은 포커스와 색감, 선명성에서 개선을 이뤄낸 것이 보인다. 최신 기종에 이를수록 개선 폭은 작아지는 것이 확인된다.
애플은 반면에 과거 기종부터 우수한 색감을 보인다. 주로 선명성과, 특히 포커싱 개선에 노력을 들인 것으로 보인다.

 

나에게 맞는 스마트폰 카메라 찾기

Moment라는 스마트폰 카메라 전문 렌즈 제작 스타트업은 “항상 휴대하는 카메라가 가장 좋은 카메라”라고 한 바 있다. 스마트폰 카메라의 우수한 휴대성을 강조한 말이다. “이왕 항상 휴대하는 거, 수치상으로 가장 좋은 카메라 대신, 내가 지향하는 사진을 찍을 수 있게 하는 카메라를 선택하는게 어떨까?”라는 생각에서 시작해서 벤치마크 숫자만으로는 파악할 수 없는 잠재적인 요인들을 요인분석으로 파악해 보았다.

색감을 중요시하는 사람은 전통적으로 우수한 색감을 자랑해온 애플 아이폰, 또는 세대가 거듭함에 따라 비약적인 색감 개선을 이뤄낸 삼성 제품을 선택할 수 있을 것이다. 또, 디테일을 중요하게 생각하는 사람은 디테일(FA1) 1위 제품을 만들어낸 구글 제품을 선택할 수 있을 것이다.

단순 벤치마크보다는 누군가의 스마트폰 선택에 도움이 되는 정보이기를 바란다.

2017년 대선후보 지지층 군집분석

대선까지 얼마 남지 않은 지금이다. 근 몇 년간 이보다 많은 사람들이 적극적으로 관심을 가진 선거가 과연 있었을까 싶다. 이 때문인지 타 후보에 대한 비방, 지지층간의 갈등과 모함 역시 이전 선거와는 궤를 달리할 정도로 적극적(?)이어 보인다. 특히나 서로 다른 후보의 지지층이 서로를 비난하는 모습은 그 정도가 지나치기도 하여, 소위 말하는 ‘국론 분열’에 대한 우려까지 자아내기도 한다.

후보자의 지지율과 그 추세, 지역, 연령대별 지지율에 대한 정보는 여론조사를 통해서 쉽게 열람할 수 있지만, 각 후보를 지지하는 지지자/지지층이 어떤 성향을 가지고 있는 지에 대해서는 쉽게 파악하기 어렵다. 그저 직관적으로 ‘이 후보를 지지하니 이런이런 사람이겠지’라고 짐작할 뿐이다. 이에 객관적으로 후보의 지지층의 성향을 파악할 수 있다면 서로에 대한 어림짐작으로부터 비롯한 비난(‘OO후보 지지하는 인간은 고령층 적폐세력 옹호자더라!’, ‘OO후보 지지자는 수도권 젊은이 뻘갱이라더라!’ 등등의 원색적인 비난 등)이 줄어들지 않을까 싶어 지지층 군집(cluster) 분석을 기획해 보았다.

 

분석 결과 이전에 알린다

  1. 후보의 성향에 대한 분석이 아니다. 후보를 지지하는 사람들의 성향에 대한 분석이다. 이 둘에는 분명한 차이가 있다.
  2. 또한 이는 후보를 지지하는 사람들의 전반적인 성향에 대한 정보를 제공하는 것이지, 지지자 개개인이 반드시 이러한 성향을 가진다는 의미를 가지지는 않는다.
  3. 모든 데이터는 중앙선거여론조사심의위원회에 결과 등록된 여론조사로, 이는 네이버 제 19대 대선 페이지에서 가져왔다. 결국 각 데이터포인트는 설문에 대한 응답 백분율(%)이다.
  4. 크게 네 가지 종류의 변수를 사용했다. 보다시피 4/7~21에 있었던 여론조사를 바탕으로 만들어진 변수이다.
    1. 정책 성향 관련
      • 일자리 정책 의견 (04.16 조선일보, 칸타퍼블릭) – 민간이 주도해야 한다고 응답한 퍼센티지
      • 사드 배치 찬반 (04.16 조선일보, 칸타퍼블릭) – 반대한다고 응답한 퍼센티지
    2. 연령대 관련
      • 세대별 ‘호감이 간다’ 응답 비율 (04.07 한국갤럽)
    3. 지역 관련
      • 지역별 후보 지지율 (04.21 동아일보, 리서치앤리서치)

    4. 지지 후보 선택 기준
      • 지지후보 결정에 중요한 기준은? (04.18 서울신문, YTN, 엠브레인) – ‘능력과 정치 경험’이라고 응답한 퍼센티지와 ‘도덕성과 청렴성’이라고 응답한 퍼센티지
      • 유승민 후보, 심상정 후보에 대한 조사가 없어서 3자 클러스터링에만 사용되었다.
  5. 데이터포인트를 클러스터링하는 데에는 계층적 클러스터링, 그 중에서도 agglomerative 방식을 사용했다. 클러스터링 옵션을 더 자세히 설명하자면 이렇다:
    • 계층적 클러스터링을 위해서는 각각의 변수값(여기서는 연령대별 지지율 등등) 사이의 유사도를 계산하는 방법인 ‘metric’과 각 샘플(즉, 여기서는 각 대선후보의 지지층) 사이의 클러스터링 방식인 ‘linkage’을 설정해주어야 한다.
    • Metric으로는 코사인 유사도(cosine similarity)를 사용했다. 이는 응답 값의 차이 대신 응답 패턴의 유사성만을 보기 위해서이다. 코사인 유사도를 사용하면 지지율과 같이 패턴이 비슷해도 값에는 차이가 클 수 있는 데이터에서 패턴의 유사도를 계산할 수 있다.
      figure 출처
    • Linkage로는 Ward’s minimum variance method를 사용했다. 그냥 많이들 사용하는 방식이라고 해서 써봤다.
    • 이 외의 metric과 linkage는 도큐멘테이션을 참고하면 좋다.
    • 클러스터링 결과로 그려진 수형도(tree)에서 가지가 짧을수록 서로 유사한 샘플(지지층)이라고 해석할 수 있다.

 

5자 지지층 클러스터링 결과

output_26_0

5자 지지층 비교 – 정책 지향

일자리 정책과 사드 배치 찬반에 대한 의견에서 정책 지향에 따라 지지층의 성향이 극명하게 나뉘는 것을 볼 수 있다. 문재인, 심상정 후보의 지지층이 정책성향이 매우 유사하고, 안철수, 유승민, 홍준표 후보의 지지층 정책성향이 매우 유사하다.

output_28_0

5자 지지층 비교 – 연령대별 호감도

연령대별 호감도로부터 주요 지지층의 연령대를 엿볼 수 있다. 20~40대과 50대 이상의 호감도가 뚜렷하게 구별됨을 확인할 수 있으며, [문재인, 심상정] 후보, [안철수, 유승민] 후보의 지지층의 연령대가 비슷할 것이라 추측할 수 있겠다. 홍준표 후보 지지층이 다른 후보 4명에 비해 유독 튀어보인다.

output_30_0

5자 지지층 비교 – 지역별 지지율

지역별로는 좌와 우 진영이 확실하게 구별되는 듯 하다. (현재는 여야랄게 없지만) 야권과 범여권의 ‘텃밭’에 따라 지지층이 확연히 구별되어 보인다.

 

항상 언론상의 뜨거운 감자인 세 후보의 지지층 정보만으로 다시 클러스터링을 해보았다. 클러스터링 옵션은 위와 같다. 다만 지지후보 선택 기준 변수가 추가되었다.

3자 지지층 클러스터링 결과

output_33_0

3자 지지층 비교 – 정책 지향

굳이 셋 만을 비교한다면 안철수 후보를 지지하는 사람들의 정책 지향은 홍준표 후보 지지자의 것과 유사도가 높게 나타난다.

output_35_0

3자 지지층 비교 – 지지후보 선택 기준

지지후보 선택 기준에 있어서는 문재인 후보 지지층과 안철수 후보 지지층이 매우 유사하다.

output_37_0

3자 지지층 비교 – 연령대별 호감도

연령대별 분류에서도 문재인 후보 지지층과 안철수 후보 지지층이 유사하다. 안철수 후보 지지층은 50대 이상에서도 꽤 나타난다는 것이 차이라면 차이.

output_39_0

3자 지지층 비교 – 지역별 지지율

지역별 분류에서도 문재인 후보 지지층과 안철수 후보 지지층이 유사하다.

 

종합하자면 이렇다

  • 정책 지향성에 있어서는 [홍준표, 유승민, 안철수] 후보군, [문재인, 심상정] 후보군으로 나눌 수 있다. 조사 대상이 두 가지 사안 뿐이었다는 한계가 있으나 사드 배치(외교, 안보)와 일자리 정책(청년, 복지, 경제)이라는 핵심 정책에 대한 조사였으므로 지지층 성향을 파악하는 데에는 충분하지 않았나 싶다.
  • 연령별 분류에 있어서는 [문재인, 심상정] / [안철수, 유승민] / [홍준표] 후보군으로 나눌 수 있다.
  • 지역별 분류에서는 [문재인, 안철수] / [홍준표, 유승민] / [심상정] 후보군으로 나눌 수 있다.
  • 문재인 후보 지지층은 전통적인 민주당계 지지층의 모습을 보여주는 듯 하다. 지역이 그러하고 연령이 그러하며 정책지향성이 그러하다. 마찬가지로 홍준표 후보 지지층은 전통적인 한나라-새누리당계 지지층의 모습을 보여주는 듯 하다.
  • 안철수 후보 지지층은 지역, 연령에 있어서는 진보로 분류되는 그룹이나, 속한 그룹과 별개로 보수적인 가치를 지지하는 모습을 보여준다. 뭇 언론에서 칭하는 것 처럼 이를 ‘중도’라 말할 수도 있겠다.

분석에 사용한 데이터 및 분석 과정과 여론조사 데이터를 긁어모으는데 사용한 스크립트는 각각 여기여기서 볼 수 있다. 실시간으로 계속해서 쌓이는 데이터여서 현재는 당시 분석에 사용한 데이터의 위치가 바뀌었을 수도 있다.

Regularized linear regression

 

import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import Lasso, Ridge, ElasticNet, LinearRegression 

X = np.random.randn(40)
y = np.cos(X) + np.random.rand(40)/3
plt.scatter(X, y)
plt.show()

output_2_0

from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
def plot_comparison(degree):
    plt.figure(figsize=(9,7))
    models = [("Linear (non-regularized)", LinearRegression()), ("Lasso", Lasso(alpha=0.01)),
              ("Ridge", Ridge(alpha=0.01)), ("Elastic Net", ElasticNet(alpha=0.01, l1_ratio=0.5))]
    for i, (model_name, model) in enumerate(models):
        plt.subplot(2,2,i+1)
        linear_model = make_pipeline(PolynomialFeatures(degree), model)
        linear_model.fit(X.reshape(-1, 1), y)
        plt.scatter(X, y, color="w", edgecolor="k")

        xx = np.linspace(-3, 4, 10000)
        plt.plot(xx, linear_model.predict(xx.reshape(-1, 1)), color="b", linewidth=1)
        plt.xlim(-4, 4)
        plt.ylim(-5, 2)
        plt.title("{} Regression".format(model_name))
    plt.tight_layout()
    plt.show()

plot_comparison(5)

output_4_0

plot_comparison(15)

output_5_0

plot_comparison(30)

output_6_0

XX = np.c_[np.random.randn(40), np.random.randn(40)**2, np.random.randn(40)**3, np.random.rand(40)*3, np.random.rand(40)]
yy = np.cos(X[0]) + np.random.rand(40)/3
alphas, coefs, dual_gaps = Lasso().path(XX, yy, alphas=np.logspace(-4, 1, 8))
plt.figure(figsize=(7,5))
for i in range(5):
    plt.plot(alphas, coefs[i], label=r"X_{}".format(i+1))
plt.axhline(0, color="k", linestyle="--")
plt.semilogx()
plt.title("Lasso path")
plt.legend()
plt.show();

output_7_0