for textmining

FastText, 실전 사용하기

06 Jul 2017 | embedding methods

이번 글에서는 페이스북에서 개발한 FastText를 실전에서 사용하는 방법에 대해 살펴보도록 하겠습니다. FastText는 구글에서 개발한 Word2Vec을 기본으로 하되 부분단어들을 임베딩하는 기법인데요. 임베딩 기법과 관련 일반적인 내용은 이곳을 참고하시면 좋을 것 같습니다.

함수 설치하기

FastText는 파이썬 gensim 패키지 내에 포함돼 주목을 받았는데요. 이상하게 제 컴퓨터 환경에서는 지속적으로 에러가 나서, 저는 페이스북에서 제공하는 C++ 기반 버전을 사용하였습니다. 이 블로그는 이 버전을 기준으로 설명할 예정입니다. 어쨌든 아래와 같은 터미널 명령어로 fastText를 내려받아 컴파일하면 바로 사용할 수 있는 상태가 됩니다.

# 설치하고자 하는 디렉토리에서 다음을 실행
$ git clone https://github.com/facebookresearch/fastText.git
$ cd fastText
$ make

이 버전은 Mac OS와 Linux 환경에서만 실행 가능합니다

입력값 준비

FastText의 입력값은 utf-8 인코딩 기반의 텍스트 파일입니다. 라인(line) 하나당 하나의 문서가 되도록 저장해 두면 됩니다. 맥을 쓰는 저는 데스크탑 폴더 밑에 fasttext 폴더를 만들어 여기에 실행파일들을 컴파일 해놓고, 데스크탑 폴더엔 다음과 같은 영화 리뷰를 ‘kor’라는 이름으로 저장해 두었습니다.

미친놈을 동경하는 이유

울었다 웃었다 종 잡을 수 없지만 끝나고 나면 2편을 찾게 되는 영화

…

어설픈북한말투 일본해표기 이제는 진짜 약빨았나 싶은 두콤비세스와 프랭코 김정은죽인것땜에그나마 별둘줌

단어벡터 학습

C++ 기반의 FastText는 터미널에서 실행 가능합니다. 학습말뭉치 파일과 결과 파일 두 가지만 지정하면 나머지는 FastText가 알아서 해줍니다. 저의 경우 맥 데스크톱 폴더에서 실행했기 때문에 맨 앞에 ‘/fastText/’라는 경로를 지정했습니다.

$ ./fastText/fasttext -input kor -output kor_model

하지만 몇 가지 추가로 옵션을 지정해주면 더욱 좋겠죠. CBOW보다는 SkipGram 모델의 성능이 나은걸로 알려져 있기 때문에 임베딩 기법은 SG를, 단어벡터의 차원수는 100을, 양옆 단어는 세개씩 보되, 말뭉치에 100번 이상 나온 단어들만 임베딩하고 싶다면 다음과 같이 실행하면 됩니다.

$ ./fastText/fasttext skipgram -input kor -output kor_model -dim 100 -ws 3 -minCount 100

다음은 FastText의 파라메터 목록입니다. 여기에서 input, output은 사용자가 지정하지 않으면 실행이 되지 않기 때문에 반드시 입력해주어야 합니다.

parameter description default
input training file path mandatory
output output file path mandatory
verbose verbosity level 2
minCount minimal number of word occurences 5
minCountLabel minimal number of label occurences 0
wordNgrams max length of word ngram 1
bucket number of buckets 2000000
minn min length of char ngram 3
maxn max length of char ngram 6
t sampling threshold 0.0001
label labels prefix []
lr learning rate 0.05
lrUpdateRate change the rate of updates for the learning rate 100
dim size of word vectors 100
ws size of the context window 5
epoch number of epochs 5
neg number of negatives sampled 5
loss loss function {ns, hs, softmax} ns
thread number of threads 12
pretrainedVectors pretrained word vectors for supervised learning []
saveOutput whether output params should be saved 0
cutoff number of words and ngrams to retain 0
retrain finetune embeddings if a cutoff is applied 0
qnorm quantizing the norm separately 0
qout quantizing the classifier 0
dsub size of each sub-vector 2

어쨌든 실행을 하면 다음과 같은 화면이 뜹니다. 학습 속도가 매우 빨라서 한국어 영화리뷰 50만여건을 학습시키는 데 맥북프로 2015 early 기준으로 5분도 채 안 걸렸던 것 같습니다.

학습결과 확인

학습이 끝나면 gensim 패키지의 Word2Vec처럼 파이썬 콘솔에서 결과를 확인할 수 있습니다. 다음과 같이 불러오면 됩니다.

from gensim.models import KeyedVectors
model = KeyedVectors.load_word2vec_format('kor')

임베딩된 단어의 리스트나 단어벡터를 뽑는 것은 다음과 같이 하면 됩니다.

# 단어 리스트 작성
vocab = model.index2word
# 전체 단어벡터 추출
wordvectors = []
for v in vocab:
	wordvectors.append(model.wv[v])

쿼리단어 기준으로 가장 유사한 단어 리스트를 뽑는 것도 gensim의 Word2Vec과 완전히 동일합니다.

# '영화'와 가장 유사한 단어 30개 뽑기
model.most_similar('영화', topn=30)

파일럿 결과

‘영화’와 가장 유사한 단어를 뽑아 봤습니다. 다음과 같습니다.

영화다, 영화였다, 작품, 영화임, 영화이다, 영화였음, 영화인듯, 영화에요, 영화네요…

‘쓰레기’와 가장 유사한 단어는 아래와 같습니다.

시간낭비, 삼류, 최악의, 최악, 아깝다, 쓰레기같은, 병신, 동, 졸작, 쓰레기를, 극혐…

Comment  Read more

베이지안 추론(2)

05 Jul 2017 | Bayes

이 글에서는 베이지안 추론에 대해 살펴보도록 하겠습니다. 이 글은 ‘Think Bayes(앨런 B. 다우니 지음, 권정민 옮김, 한빛미디어 펴냄)’를 정리했음을 먼저 밝힙니다. 그럼 시작하겠습니다.

기관차 문제

모든 기차의 기관차엔 1부터 $N$까지 일련번호가 붙어 있습니다. 어느 날 60호 기관차를 봤다고 칩시다. 이 때 이 회사가 운영하고 있는 기관차는 총 몇 개일까요? 즉 $N$이 몇인지 추정하는 것이 문제가 되겠습니다. 베이지안 추론법에 따르면 이런 문제는 두 단계로 해결합니다.

(1) 사전확률 : 데이터를 보기 전에 $N$에 대해 알고 있는 것이 무엇인가?

(2) 우도 : $N$이 특정한 값을 가질 때 관측한 데이터(60호 객차)가 발생할 확률은?

동일 사전확률 가정

사전확률과 관련해 근거는 부족하지만 간단한 가정을 해보겠습니다. 철도를 운영하는 회사는 단 하나이고, 이 회사가 운영 중인 기관차는 일단 1000개라고 칩시다. 또 $N$은 1~1000 사이의 어떤 값이든 동일한 확률로 선택될 수 있다고 보는 겁니다.

이 가정 하에서 우도를 구해봅시다. 이 회사가 운영하고 있는 기관차가 59개 이하라면 60호 기관차를 절대 관측할 수 없습니다. 따라서 $N$이 1~59일 때 우도는 0입니다.

운영 중인 기관차 수가 60개라면 60호 기관차를 관측할 확률은 1/60이 됩니다. 마찬가지로 총 기관차 수가 61개라면 우도는 1/61, 62개라면 1/62… 이런 식으로 우도를 구할 수 있게 됩니다.

베이즈 규칙에 따르면 사후확률(60호 기관차를 관측했을 때 운영 중인 기관차 수 $N$이 나타날 확률)은 사전확률에 우도를 곱해 풀 수 있습니다. 사전확률은 동일하다고 가정했기 때문에 제일 큰 값을 갖는 우도가 전제하는 $N$이 우리의 추정 결과가 될 것입니다.

이 경우 우리는 $N$을 60으로 추정할 수밖에 없게 됩니다. 운영 중인 기관차가 60개일 때 우도가 가장 큰 값을 갖게 되거든요. 그런데 이런 추정은 상식에는 반하는 결과입니다.

다시 말해 60호 기관차를 관측했다고 해서 이 회사가 운영 중인 기관차 수가 60개라고 추정하는 건 너무 단순한 추론이라는 것이지요. 만약 35호 기관차를 관측했다면 같은 방식으로 베이지안 추론을 해보면 총 기관차 수가 35개라는 예측이 나오게 될 겁니다.

이 가정에 따른 결과 표는 다음과 같습니다. 아래 표에서 사후확률은 사전확률에 우도를 곱한 뒤 사후확률 전체 합이 1이 되도록 나누어준 값입니다.

$N$ 사전확률 우도 사후확률
1 1/1000 0 0
… … …  
59 1/1000 0 0
60 1/1000 1/60 0.0059
… … …  
1000 1/1000 1/1000 0.00035

사후확률의 평균은 $N$에 사후확률을 곱한 뒤 모두 더해주면 됩니다. 그 값은 333.065557입니다.

가정 보완하기

동일 사전확률 가정이 불완전함을 확인했습니다. 이를 보완할 수 있는 방법은 두 가지가 있습니다. (1)이 사전확률, (2)가 우도 쪽을 보완하는 접근기법이 되겠습니다.

(1) 배경지식을 더 확보

(2) 데이터를 더 확보

사전확률의 대안

(1)과 관련해 기관차를 1000대나 운영하고 있는 회사가 있고, $N$이 1~1000 사이에 선택될 확률이 동일하다고 전제하는 건 비현실적입니다. 회사의 규모는 현실적으로는 이보다 규모가 훨씬 작을 겁니다. 바꿔 말해 $N$이 1000이 될 확률은 $N$이 1이 될 가능성보다 훨씬 작을 것이라 가정하는 것이 합리적입니다.

이와 관련해 로버트 액스텔이라는 학자는 회사 규모의 분포는 멱법칙을 따른다는 사실을 밝혀냈습니다. 이 법칙에 따르면 기관차가 10대 미만인 회사가 1000개일 때, 100대의 기관차를 소유한 회사는 100개, 1000대의 기관차를 운영 중인 회사는 10개, 1만대 기관차를 소유한 회사는 1개입니다.

다시 말해 운영 중인 기관차의 총 개수 $N$이 커질수록 $N$이 선택될 확률, 즉 사전확률이 점차 작아진다는 이야기입니다. 그럼 사전확률을 다시 구해보겠습니다.

우선 $N$의 역수를 취합니다.

$N$ 역수
1 1
2 0.5
3 0.33
4 0.25
5 0.2
… …
총합 7.48547

이렇게 구한 역수를 역수의 총합(7.48547)으로 나눠준 것이 각 $N$의 사전확률이 되겠습니다. 다음과 같습니다.

$N$ 사전확률
1 0.13359213049244018
2 0.06679606524622009
3 0.044530710164146725
4 0.033398032623110044
5 0.026718426098488037
… …
996 0.0001341286450727311
997 0.00013399411283093298
998 0.00013385985019282582
999 0.00013372585634878898
1000 0.00013359213049244018

관측치 늘리기

동일 사전확률 가정 하에서 60호 기관차 하나만 관측한 뒤에 사후확률의 평균을 구하면 다음과 같습니다.

상한선 사후확률 평균
500 207
1000 333
2000 552

우리는 60호 기관차만 보고 결론을 내린 셈인데요. 데이터를 하나만 보고 결론을 내리는건 너무 성급한 일입니다. 60호 기관차에 이어 30번과 90번 기관차도 봤다고 가정해 보겠습니다. 이 데이터를 사용했을 때 사후확률 평균은 다음과 같습니다.

상한선 사후확률 평균
500 152
1000 164
2000 171

그 차이가 확실히 줄어들어 수렴하는 경향을 보임을 알 수 있습니다.

마지막으로 보완한 사전확률을 토대로 30번, 60번, 90번 기관차를 관측한 뒤의 사후확률 평균을 보겠습니다.

상한선 사후확률 평균
500 131
1000 133
2000 134

상한선이 어떻든 간에 상관없이 사후확률 평균은 134에 수렴한다고 합니다.

파이썬 코드

마지막 사후확률 표를 생성하는 데 쓴 파이썬 코드는 다음과 같습니다. thinkbayes.py는 이곳에서 내려받을 수 있습니다.

import thinkbayes as tb

class Train(tb.Suite):
    """Represents hypotheses about 
    how many trains the company has."""

    def __init__(self, hypos, alpha=1.0):
        """Initializes the hypotheses 
        with a power law distribution.

        hypos: sequence of hypotheses
        alpha: parameter of the power law prior
        """
        tb.Pmf.__init__(self)
        for hypo in hypos:
            self.Set(hypo, hypo ** (-alpha))
        self.Normalize()

    def Likelihood(self, data, hypo):
        """Computes the likelihood of the data 
        under the hypothesis.
        """
        if hypo < data:
            return 0
        else:
            return 1.0 / hypo

def Mean(suite):
    total = 0
    for hypo, prob in suite.Items():
        total += hypo * prob
    return total

def MakePosterior(high, dataset):
    hypos = xrange(1, high+1)
    suite = Train(hypos)
    suite.name = str(high)

    for data in dataset:
        suite.Update(data)

    return suite

dataset = [30, 60, 90]
for high in [500, 1000, 2000]:
    suite = MakePosterior(high, dataset)
    print high, suite.Mean()

Comment  Read more

선형회귀 파라메터 추정

03 Jul 2017 | regression

이번 글에서는 선형회귀 모델의 계수를 추정하는 방법을 살펴보도록 하겠습니다. 이번 글은 고려대 김성범 교수님 강의와 ‘밑바닥부터 시작하는 데이터과학(조엘 그루스 지음, 인사이트 펴냄)’을 정리하였음을 먼저 밝힙니다. 그럼 시작하겠습니다.

선형회귀

선형회귀(Multiple Linear Regression)는 수치형 설명변수 X와 연속형 숫자로 이뤄진 종속변수 Y간의 관계를 선형으로 가정하고 이를 가장 잘 표현할 수 있는 회귀계수를 데이터로부터 추정하는 모델입니다. 다음 그림처럼 집 크기와 가격의 관계를 나타내는 직선을 찾는 것이 선형회귀의 목표입니다.

독립변수들로 이뤄진 행렬 $X$와 종속변수 벡터 $Y$가 주어졌을 때 다중선형회귀 모델은 다음과 같이 정의됩니다.

[\begin{bmatrix} { y }{ 1 } \ { y }{ 2 } \ … \ { y }{ n } \end{bmatrix}=\begin{bmatrix} 1 & x{ 11 } & … & { x }{ 1k } \ 1 & { x }{ 21 } & … & { x }{ 2k } \ … & … & & … \ 1 & { x }{ n1 } & … & { x }{ nk } \end{bmatrix}\begin{bmatrix} { \beta }{ 0 } \ { \beta }{ 1 } \ … \ { \beta }{ k } \end{bmatrix}+\begin{bmatrix} { \varepsilon }{ 1 } \ { \varepsilon }{ 2 } \ … \ { \varepsilon }_{ n } \end{bmatrix}\ \ \overrightarrow { y } =X\overrightarrow { \beta } +\overrightarrow { \varepsilon } ,\quad \overrightarrow { \varepsilon } \sim N(E(\overrightarrow { \varepsilon } ),V(\overrightarrow { \varepsilon } ))\ \ E(\overrightarrow { \varepsilon } )=\begin{bmatrix} 0 \ 0 \ … \ 0 \end{bmatrix},\quad V(\overrightarrow { \varepsilon } )={ \sigma }^{ 2 }I]

Direct Solution

선형회귀의 계수들은 실제값과 모델 예측값의 차이, 즉 오차제곱합(error sum of squares)을 최소로 하는 값들입니다. 이를 만족하는 최적의 계수들은 회귀계수에 대해 미분한 식을 0으로 놓고 풀면 아래와 같이 명시적인 해를 구할 수 있습니다. 다시 말해 우리에게 주어진 $X$, $Y$ 데이터만 가지고 계수를 단번에 추정할 수 있다는 이야기입니다.

[\overrightarrow { \beta } ={ \left( { X }^{ T }X \right) }^{ -1 }{ X }^{ T }\overrightarrow { y }]

분석 대상 데이터

분석 대상 데이터는 ‘밑바닥부터 시작하는 데이터 과학’에서 제시된 예시데이터입니다. $X$는 친구수, 근무시간, 박사학위 취득여부이고, $Y$는 사용자가 웹사이트에서 보내는 시간(분)에 해당합니다. $X$ 가운데 친구 수 변수 하나만 떼어 $Y$와의 관계를 2차원으로 도시하면 아래 그림과 같습니다.

다음은 데이터입니다. $X$는 네 개 요소로 구성돼 있는데 각각 상수항, 친구수, 근무시간, 박사학위 취득 여부에 해당합니다.

# [1(beta_0), 친구수, 근무시간, 박사학위 취득 여부]
input = [[1,49,4,0],[1,41,9,0],[1,40,8,0],[1,25,6,0],[1,21,1,0],[1,21,0,0],[1,19,3,0],[1,19,0,0],[1,18,9,0],[1,18,8,0],[1,16,4,0],[1,15,3,0],[1,15,0,0],[1,15,2,0],[1,15,7,0],[1,14,0,0],[1,14,1,0],[1,13,1,0],[1,13,7,0],[1,13,4,0],[1,13,2,0],[1,12,5,0],[1,12,0,0],[1,11,9,0],[1,10,9,0],[1,10,1,0],[1,10,1,0],[1,10,7,0],[1,10,9,0],[1,10,1,0],[1,10,6,0],[1,10,6,0],[1,10,8,0],[1,10,10,0],[1,10,6,0],[1,10,0,0],[1,10,5,0],[1,10,3,0],[1,10,4,0],[1,9,9,0],[1,9,9,0],[1,9,0,0],[1,9,0,0],[1,9,6,0],[1,9,10,0],[1,9,8,0],[1,9,5,0],[1,9,2,0],[1,9,9,0],[1,9,10,0],[1,9,7,0],[1,9,2,0],[1,9,0,0],[1,9,4,0],[1,9,6,0],[1,9,4,0],[1,9,7,0],[1,8,3,0],[1,8,2,0],[1,8,4,0],[1,8,9,0],[1,8,2,0],[1,8,3,0],[1,8,5,0],[1,8,8,0],[1,8,0,0],[1,8,9,0],[1,8,10,0],[1,8,5,0],[1,8,5,0],[1,7,5,0],[1,7,5,0],[1,7,0,0],[1,7,2,0],[1,7,8,0],[1,7,10,0],[1,7,5,0],[1,7,3,0],[1,7,3,0],[1,7,6,0],[1,7,7,0],[1,7,7,0],[1,7,9,0],[1,7,3,0],[1,7,8,0],[1,6,4,0],[1,6,6,0],[1,6,4,0],[1,6,9,0],[1,6,0,0],[1,6,1,0],[1,6,4,0],[1,6,1,0],[1,6,0,0],[1,6,7,0],[1,6,0,0],[1,6,8,0],[1,6,4,0],[1,6,2,1],[1,6,1,1],[1,6,3,1],[1,6,6,1],[1,6,4,1],[1,6,4,1],[1,6,1,1],[1,6,3,1],[1,6,4,1],[1,5,1,1],[1,5,9,1],[1,5,4,1],[1,5,6,1],[1,5,4,1],[1,5,4,1],[1,5,10,1],[1,5,5,1],[1,5,2,1],[1,5,4,1],[1,5,4,1],[1,5,9,1],[1,5,3,1],[1,5,10,1],[1,5,2,1],[1,5,2,1],[1,5,9,1],[1,4,8,1],[1,4,6,1],[1,4,0,1],[1,4,10,1],[1,4,5,1],[1,4,10,1],[1,4,9,1],[1,4,1,1],[1,4,4,1],[1,4,4,1],[1,4,0,1],[1,4,3,1],[1,4,1,1],[1,4,3,1],[1,4,2,1],[1,4,4,1],[1,4,4,1],[1,4,8,1],[1,4,2,1],[1,4,4,1],[1,3,2,1],[1,3,6,1],[1,3,4,1],[1,3,7,1],[1,3,4,1],[1,3,1,1],[1,3,10,1],[1,3,3,1],[1,3,4,1],[1,3,7,1],[1,3,5,1],[1,3,6,1],[1,3,1,1],[1,3,6,1],[1,3,10,1],[1,3,2,1],[1,3,4,1],[1,3,2,1],[1,3,1,1],[1,3,5,1],[1,2,4,1],[1,2,2,1],[1,2,8,1],[1,2,3,1],[1,2,1,1],[1,2,9,1],[1,2,10,1],[1,2,9,1],[1,2,4,1],[1,2,5,1],[1,2,0,1],[1,2,9,1],[1,2,9,1],[1,2,0,1],[1,2,1,1],[1,2,1,1],[1,2,4,1],[1,1,0,1],[1,1,2,1],[1,1,2,1],[1,1,5,1],[1,1,3,1],[1,1,10,1],[1,1,6,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],[1,1,4,1],[1,1,9,1],[1,1,9,1],[1,1,4,1],[1,1,2,1],[1,1,9,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],[1,1,1,1],[1,1,1,1],[1,1,5,1]]
# 사이트에서 보내는 시간(분)
output = [68.77,51.25,52.08,38.36,44.54,57.13,51.4,41.42,31.22,34.76,54.01,38.79,47.59,49.1,27.66,41.03,36.73,48.65,28.12,46.62,35.57,32.98,35,26.07,23.77,39.73,40.57,31.65,31.21,36.32,20.45,21.93,26.02,27.34,23.49,46.94,30.5,33.8,24.23,21.4,27.94,32.24,40.57,25.07,19.42,22.39,18.42,46.96,23.72,26.41,26.97,36.76,40.32,35.02,29.47,30.2,31,38.11,38.18,36.31,21.03,30.86,36.07,28.66,29.08,37.28,15.28,24.17,22.31,30.17,25.53,19.85,35.37,44.6,17.23,13.47,26.33,35.02,32.09,24.81,19.33,28.77,24.26,31.98,25.73,24.86,16.28,34.51,15.23,39.72,40.8,26.06,35.76,34.76,16.13,44.04,18.03,19.65,32.62,35.59,39.43,14.18,35.24,40.13,41.82,35.45,36.07,43.67,24.61,20.9,21.9,18.79,27.61,27.21,26.61,29.77,20.59,27.53,13.82,33.2,25,33.1,36.65,18.63,14.87,22.2,36.81,25.53,24.62,26.25,18.21,28.08,19.42,29.79,32.8,35.99,28.32,27.79,35.88,29.06,36.28,14.1,36.63,37.49,26.9,18.58,38.48,24.48,18.95,33.55,14.24,29.04,32.51,25.63,22.22,19,32.73,15.16,13.9,27.2,32.01,29.27,33,13.74,20.42,27.32,18.23,35.35,28.48,9.08,24.62,20.12,35.26,19.92,31.02,16.49,12.16,30.7,31.22,34.65,13.13,27.51,33.2,31.57,14.1,33.42,17.44,10.12,24.42,9.82,23.39,30.93,15.03,21.67,31.09,33.29,22.61,26.89,23.48,8.38,27.81,32.35,23.84]

$i$번째 데이터 포인트에 대한 회귀식은 다음과 같습니다. 우리는 여기에서 회귀계수로 이뤄진 벡터 $β$를 추정해야 합니다.

[{ y }{ i }=\beta { x }{ i }]

$(X^TX)^{-1}X^Ty$를 계산해 데이터로부터 명시적인 해를 구한 결과 $β$는 [30.63, 0.97, -1.87, 0.91]로 추정되었습니다.

경사하강법 같은 반복적인 방식으로 선형회귀 계수를 구할 수도 있습니다. 경사하강법이란 어떤 함수값을 최소화하기 위해 임의의 시작점을 잡은 후 해당 지점에서의 그래디언트(경사)를 구하고, 그래디언트의 반대 방향으로 조금씩 이동하는 과정을 여러번 반복하는 것입니다. 예컨대 아래 그림(출처)과 같습니다.

이 글에서는 경사하강법 가운데 Stochastic Gradient Descent(SGD) 기법을 쓰겠습니다. SGD는 반복문을 돌 때마다 개별 데이터 포인트에 대한 그래디언트를 계산하고 이 그래디언트의 반대 방향으로 파라메터를 업데이트해 함수의 최소값을 구하는 기법입니다.

SGD 관련 메인 코드는 다음과 같습니다. 코드에서 ‘value’는 우리가 최소화하고 싶은 값으로 선형회귀 모델에서는 오차제곱합을 가리킵니다. ‘target_fn’은 목적함수로 오차제곱합을 아웃풋으로 산출하는 함수를 의미합니다. ‘theta’는 해당 목적함수의 파라메터인데요, 우리 문제에선 $α$와 $β$를 말합니다. ‘gradient_fn’은 각 파라메터에 대한 목적함수의 그래디언트를 가리킵니다.

def minimize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
    # SGD 방식으로 gradient descent
    # minimize_batch보다 훨씬 빠르다

    data = zip(x, y)
    # theta_0를 초기 theta로
    theta = theta_0
    # alpha_0를 초기 이동거리(step_size)로
    alpha = alpha_0
    # 시작할 때의 최소값
    min_theta, min_value = None, float("inf")
    iterations_with_no_improvement = 0

    # 만약 100번 넘게 반복하는 동안 value가 더 작아지지 않으면 멈춤
    while iterations_with_no_improvement < 100:
        value = sum( target_fn(x_i, y_i, theta) for x_i, y_i in data )

        # 새로운 최솟값을 찾았다면
        if value < min_value:
            # 이 값을 저장
            min_theta, min_value = theta, value
            # 100번 카운트도 초기화
            iterations_with_no_improvement = 0
            # 기본 이동거리로 돌아감
            alpha = alpha_0

        # 만약 최솟값이 줄어들지 않는다면
        else:
            # 이동거리 축소
            alpha *= 0.9
            # 100번 카운트에 1을 더함
            iterations_with_no_improvement += 1
		
        # 반복문이 돌 때마다 in_random_order를 호출하기 때문에
        # 매 iter마다 그래디언트를 계산하는 순서가 달라짐
        for x_i, y_i in in_random_order(data):
            # 각 데이터 포인트에 대해 그래디언트를 계산
            gradient_i = gradient_fn(x_i, y_i, theta)
            # 기존 theta에서, 학습률(alpha)과 그래디언트를 뺀 것을 업데이트
            theta = vector_subtract(theta, scalar_multiply(alpha, gradient_i))

    return min_theta

그러면 이제는 value와 target_fn, gradient_fn, theta를 정의해야 합니다. 우선 우리가 구해야 하는 파라메터는 입력변수 개수($k$=3)+상수항(1) 길이를 가진 벡터 $β$이므로 beta는 4차원 벡터로 선언했습니다.

import random
random.seed(10)
# 추정 대상 beta (vector)
# beta = [beta_0, beta_1, ..., beta_k]
beta = [random.random() for x_i in input[0]]

우리가 최소화하고자 하는 값(value)은 오차제곱합입니다. $i$번째 데이터 포인트에 대한 오차제곱(Squared Error)은 다음과 같은 식으로 나타낼 수 있습니다.

[{SE}{i}={ \left{ { y }{ i }-\left( \beta { x }_{ i }\right) \right} }^{ 2 }]

이를 ‘squared_error’ 함수로 표현할 수 있습니다.

def squared_error(x_i, y_i, beta):
    return error(beta, x_i, y_i) ** 2

def error(beta, x_i, y_i):
    # 실제값 y_i와 예측값 사이의 편차
    return y_i - predict(beta, x_i)

def dot(v, w):
    """v_1 * w_1 + ... + v_n * w_n"""
    return sum(v_i * w_i for v_i, w_i in zip(v, w))

def predict(beta, x_i):
    # 현재 회귀계수들을 가지고 예측
    # 예측값 = x_i와 beta의 선형결합
    # x_i = [1(beta_0), x_i1, x_i2, ..., x_ik]
    # beta = [beta_0, beta_1, ..., beta_k]
    return dot(x_i, beta)

gradient_fn은 다음과 같습니다. 목적함수인 ‘squared_error’를 ‘theta’로 미분한 값입니다. 이를 식으로 정리하면 다음과 같습니다.

[\frac { \partial { SE }{ i } }{ \partial \alpha } =-2\left{ { y }{ i }-\left( \beta { x }{ i }+\alpha \right) \right} \ \frac { \partial { SE }{ i } }{ \partial \beta } =-2\left{ { y }{ i }-\left( \beta { x }{ i }+\alpha \right) \right} { x }_{ i }]

이를 코드로 나타내면 다음과 같습니다.

def squared_error_gradient(x_i, y_i, beta):
    # i번째 오류제곱 값의 beta에 대한 편미분값
    return [-2 * x_ij * error(beta, x_i, y_i) for x_ij in x_i]

이밖에 ‘minimize_stochastic’ 구동에 필요한 함수도 정의하겠습니다.

def in_random_order(data):
    # Stochastic Gradient Descent 수행을 위한 함수로,
    # 한번 반복문을 돌 때마다 임의의 순서로 데이터 포인트를 반환
    # 데이터 포인트의 인덱스를 list로 생성
    indexes = [i for i, _ in enumerate(data)]
    # 이 인덱스를 랜덤하게 섞는다
    random.shuffle(indexes)
    # 이 순서대로 데이터를 반환한다
    for i in indexes:
        yield data[i]

def vector_subtract(v, w):
    """subtracts two vectors componentwise"""
    return [v_i - w_i for v_i, w_i in zip(v,w)]

def scalar_multiply(c, v):
    return [c * v_i for v_i in v]

마지막으로 코드 전체를 구동하는 명령문은 다음과 같습니다.

beta = minimize_stochastic(target_fn=squared_error,
                           gradient_fn=squared_error_gradient,
                           x=input,
                           y=output,
                           theta_0=beta,
                           alpha_0=0.0001)

SGD 기법으로 해를 구한 결과 $β$는 [30.55184071133236, 0.973395629801253, -1.8632386392995979, 0.9471533590851985]로 추정되었습니다. 이는 명시적 해와 유사합니다.

회귀계수를 얼마나 신뢰할 수 있나

데이터로부터 추정한 회귀계수는 진짜 회귀계수라고 할 수는 없습니다. 어디까지나 일부 데이터를 가지고 도출된 계수일 뿐더러 데이터에 노이즈가 끼어있을 수도 있기 때문이죠. 선형회귀 모델의 $j$번째 독립변수에 대한 추정 회귀계수 $β_j$를 추정 회귀계수의 표준편차 $σ_j$로 나눈 값($t_j$)은 $n-k$의 자유도를 지닌 $t$분포를 따른다는 사실이 알려져 있습니다.

그렇다면 회귀계수의 표준편차는 어떻게 구할까요? 이 때 쓰는 것이 bootstap입니다. 기존 데이터에서 중복을 허용된 재추출을 통해 새로운 데이터를 만들어내는 방법입니다. 우리가 갖고 있는 학습데이터에 bootstrap을 적용하여 여러 새로운 데이터를 생성하고, 이들 각각으로부터 추정 회귀계수를 도출할 수 있습니다.

예를 들어 특정 독립변수에 해당하는 회귀계수가 부트스트랩 데이터마다 크게 달라지지 않는다면 해당 회귀계수는 상당히 신뢰할 수 있을 겁니다. 반대로 계수가 크게 변한다면 추정된 계수는 신뢰할 수 없을 겁니다.

부트스트랩 데이터 개수만큼 회귀계수 추정값을 구하여 회귀계수의 표준편차 $σ_j$를 구한 뒤 이를 원래 데이터에서 구한 추정회귀 계수 $β_j$에 나눠줘 $t_j$를 계산합니다. 이 $t_j$가 $t$분포 상에서 어느 위치에 해당하는지를 따져서 진짜 회귀계수가 속해 있을 수 있는 신뢰구간을 구할 수 있게 됩니다.

부트스트랩 데이터 각각에 대해 회귀계수를 먼저 구해보겠습니다. 우선 데이터로부터 회귀계수를 추정하는 함수를 만들었습니다.

def estimate_beta(x, y):
    # 추정 대상 beta (vector)
    beta = [random.random() for _ in x[0]] 
    return minimize_stochastic(target_fn=squared_error,
                           gradient_fn=squared_error_gradient,
                               x=x,
                               y=y,
                               theta_0=beta,
                               alpha_0=0.0001)

이제는 부트스트랩 함수를 만들 차례입니다.

def bootstrap_sample(data):
    # 데이터 개수만큼 data로부터 무작위 재추출 (중복 허용)
    return [random.choice(data) for _ in data]

def bootstrap_statistic(data, stats_fn, num_samples):
    # num_sample개의 bootstrap 샘플에 대해 stats_fn을 적용
    return [stats_fn(bootstrap_sample(data))
            for _ in range(num_samples)]

def estimate_sample_beta(sample):
    # sample은 (x_i, y_i)로 구성된 리스트
    x_sample, y_sample = zip(*sample)
    return estimate_beta(x_sample, y_sample)

이제 함수를 실행해 봅시다.

bootstrap_betas = bootstrap_statistic(zip(input, output),
                                      estimate_sample_beta,
                                      100)

함수를 실행하면 100개의 부트스트랩 데이터에 대해 추정된 회귀계수의 쌍이 100개가 도출됩니다. 다음과 같습니다.

[[30.343293717263954, 0.955828442766452, -1.8249710660990814, 0.19262090720112343], 
[31.653090581087664, 0.8930047770677926, -1.8274769383940832, 0.003915653745906352], 
[30.062130836429326, 0.9735595051657044, -1.7409587831575586, -0.00993645069080973], 
[30.459536927917853, 0.9275688208209533, -1.970156796143431, 1.4162870885210503], 
[29.965245738566743, 0.9987255912413002, -1.8469983707870146, 1.2570210865880311], 
[30.688436939427994, 0.949721374789511, -1.791565743167013, 0.8490708882537681], 
[27.389069317384713, 1.1899927371032066, -1.6395966566222984, 2.201209799844561], 
[31.578900415746247, 0.963470036264341, -2.061237515990318, 0.6044582651703277], 
[29.9969538479282, 1.0563305743635187, -2.0370138317358597, 1.395128459047288], 
[30.71129774327976, 0.9412707718288738, -1.9049341856989015, 0.9975476331019565], 
[28.95704337705777, 1.0372235777797778, -1.7374326711464978, 2.017330974451086], 
[29.347442580389554, 0.9691244573476897, -1.6980453130961777, 1.758299750670592], 
[30.89303821117151, 0.9643400467987878, -2.0014985101424188, 0.9437323561405635], 
[31.033963018013235, 0.944512876776604, -1.8339869584884239, 1.3762002169288328], 
[29.7724071122711, 0.9637788276539256, -1.8712203112136903, 2.31343127588546], 
[30.712161312804817, 0.9787589052763698, -1.9637325222425992, 1.1877691857117059], 
[28.882828505113785, 1.0615985715436485, -1.8322410460473135, 2.282737207651317], 
[31.164357855104715, 0.9313158660935656, -1.9362632274947122, 1.190519238545326], 
[31.02991081861022, 0.9646849000104669, -2.048540894739888, 1.4641504377576209], 
[31.270216131821417, 0.9412983344101367, -1.9754923012624737, 0.9244505522676154], 
[28.972811257580727, 1.0626659243331988, -1.8172782368438227, 1.9772010790265988], 
[32.47326463848947, 0.9446160353528085, -1.8827511459492636, -0.555442587406212], 
[31.92477575189198, 0.9076629023019004, -1.833715795651844, 0.2543582162753717], 
[30.73653173338487, 0.9798526580610298, -1.8368633847870894, 1.3991627628282644], 
[31.052387396870795, 0.9465884318410379, -1.937853405314178, 1.6167719138597119], 
[30.838191991050294, 0.9547770115578801, -1.9540199662139435, 1.1143266882659546], 
[30.989458366759308, 0.9112712446387848, -1.8275575068998577, 1.1518607341417386], 
[31.499881551689803, 0.984790498881108, -2.017772833184302, 0.6959781050813002], 
[30.6294509899701, 0.9532711074826222, -1.9002085912683446, 0.41904632258007973], 
[31.871492871773135, 0.9640564217099683, -2.0592781960142137, -0.1325017160959066], 
[30.979163018183787, 0.9442339015956216, -1.9852723070949032, 1.018733620546921], 
[30.94960897962424, 0.9917550347995638, -1.879826427606344, -0.9344395120511418], 
[29.974194993513784, 0.9372566526240828, -1.6983466191579584, 1.413895309854322], 
[30.48871554272555, 0.9274603561875837, -1.6730463683580687, 0.7292989104575046], 
[29.976667988625856, 0.9711827384011987, -1.839308476190981, 1.9172708899219006], 
[30.692522310161685, 0.951972668383523, -1.7681190534968108, 0.4691416207459495], 
[28.288129698240517, 1.0948205153184756, -1.8171117444664553, 3.3038823592912854], 
[29.822963056351064, 1.0327653743440959, -1.869836637120832, 1.5486447223042832], 
[30.386808122240165, 0.9779154202444837, -1.8976774910373846, 2.1086881935517914], 
[30.741619085348518, 0.9217376178433008, -1.7198893947497884, -0.34894837617535585], 
[29.12768874650376, 1.0267526931238684, -1.7647591340646893, 1.6875260877074156], 
[30.24960760356648, 1.0121231682791874, -1.7376567600056418, 0.296763706300862], 
[31.035301020344654, 0.9431366335375586, -1.953384748218975, 1.1579872338326034], 
[30.5768817821496, 0.947021358453022, -1.7205109807167238, 0.783267985369746], 
[29.374898109259238, 1.0939268639920088, -1.6424719680107713, 1.140410049137131], 
[29.40635067130438, 1.0637366346977501, -1.874508187809908, 2.4511133687480937], 
[30.95063057843148, 0.936325594647627, -1.7737167596585748, 1.0343706367628467], 
[31.678274244399336, 0.9227541184400466, -1.9918345207726629, 0.5086225151581116], 
[29.98166257475079, 0.9605170287193446, -1.7577319419624025, 0.5120422638922413], 
[31.323616890901306, 0.8176945107731992, -1.6500326316218685, -0.8683616760822116], 
[30.302149825423008, 0.9737686518258278, -1.8128215160225314, 0.7033976674704967], 
[31.386502073193892, 0.882658769388189, -1.6825286945756024, -0.44420246297894167], 
[30.323858845499185, 1.0757705725166067, -2.071948018422145, 2.2968801639230105], 
[30.494825872775092, 1.0108387224084352, -1.9188457540102233, 0.8717781457244161], 
[30.576615201544183, 0.9120017225776094, -1.7251677135358494, 1.7832767592542957], 
[28.492384747280074, 1.1059593661729585, -1.7735172156402836, 2.663999434050509], 
[30.104972970041402, 1.0840732600063667, -1.9613608117473043, 1.5315485478845854], 
[31.757239201637674, 0.9596734272864389, -2.051685205450117, -0.2098820060083817], 
[29.854863672046072, 1.0209262822106175, -1.9439212461621536, 1.5185447249181654], 
[29.41528261809298, 1.1023172226932805, -1.891600842305552, 2.2231352104855175], 
[29.11300285931626, 1.0343057536723956, -1.7718637081275448, 1.6750740838338922], 
[29.39039719035395, 0.9297526750567007, -1.7896140395445401, 1.940908892848941], 
[27.962469288747435, 1.1422689401264716, -1.8893543588108115, 3.388049453882638], 
[29.099157058356646, 1.1348583693290986, -1.8317976244428604, 1.663037801381304], 
[29.137612379810427, 1.1300566948876725, -1.7913898643267325, 1.6482641446500326], 
[31.130222918185655, 0.9030219758748801, -1.9131692904927817, 1.474759950650903], 
[30.58125248343542, 0.9223235418040419, -1.6101280249243748, 1.0127434495055474], 
[30.891973319804126, 0.9193557426060519, -1.8375241356871626, 0.08461128745246431], 
[31.14771990918551, 0.9553305770320748, -1.878558304343514, 0.9002047428943074], 
[31.27678282839363, 0.9335965433539816, -1.9467073427528274, 0.9226402524279342], 
[29.56185170659549, 1.077105492914354, -1.8241298785196545, 1.595589221533047], 
[28.076434296423237, 1.0481641153277417, -1.6920455897818678, 1.8737381674617253], 
[28.90380434686726, 1.0604153971116186, -1.7365915272861396, 2.026884447765734], 
[29.99741704319044, 1.0666687600651432, -1.9111028667834695, 1.2256669988576885], 
[29.43736242616701, 1.0472064322799974, -1.8063134614295004, 1.9580133762596552], 
[29.227650225513607, 1.0144255050958706, -1.8016989061013566, 1.7520195371651819], 
[29.243279432899033, 0.9968822418447205, -1.778030031899557, 0.5551593556261414], 
[31.46888932259628, 0.9105043256707787, -1.7449900371478968, -0.7085272038815223], 
[29.924868298448885, 1.0239101298478512, -1.8708758584829717, 0.959971954552992], 
[30.88467060799466, 0.9380489354101933, -1.8289315911516493, 0.36156696066090643], 
[31.757648253366902, 0.96456365074318, -2.12281082826505, 0.4933404740397043], 
[30.359636241287728, 0.9930935207021235, -1.91614726756068, 0.11981014378409671], 
[29.689503027235425, 1.0275941570260614, -1.8684278630889204, 1.695936319407561], 
[30.675188326985708, 1.045569865938166, -1.9146598992756043, 1.0011598992584545], 
[30.43186544380135, 1.0040689827077782, -2.0168868872700227, 1.3592129670325508], 
[31.56317410958609, 0.9256366346831482, -1.8198369687354612, 0.37814324357003987], 
[32.0724948565735, 0.9134883775996528, -1.9211327000011957, -1.1653914277078095], 
[29.24454398282765, 0.9972606107989663, -1.7192410689836413, 1.0773787739745337], 
[29.565133763140803, 0.969925575012797, -1.6784388902600302, 1.16934425282509], 
[30.521873388440795, 0.9659927881039246, -1.81299769918763, 0.9066914988998109], 
[31.304315710353873, 0.9388059361772737, -1.7375886505932063, 0.22777427755049626], 
[30.485899853898623, 0.9831955516685744, -1.8216915684007813, 1.4665292532815406], 
[30.586128743092342, 0.954830980894403, -1.9606989544282951, 2.7599642168405083], 
[30.331607238578425, 0.9800674077114176, -1.9621069206081838, 2.229119033504103], 
[31.372326222839966, 1.0047554827884049, -2.0193398364534048, 0.9477139083724562], 
[30.231048927885226, 1.0141246515657603, -1.8918513233859755, 0.1924866352716917], 
[30.263326574321635, 0.9311507918160458, -1.9097048979881397, 2.9233469946337465], 
[31.515503088339216, 0.9178784037775057, -1.998991069549314, 1.2935446029780655], 
[30.03038991481507, 1.035274803835536, -2.028483027984759, 2.544621079119659], 
[31.070081646523697, 0.927782273878715, -1.9091320174841293, 0.15186481410364247]]

우리는 이로부터 회귀계수 각각의 표준편차를 구할 수 있습니다. SGD 기법으로 해를 구한 결과 $β$는 [30.55184071133236, 0.973395629801253, -1.8632386392995979, 0.9471533590851985]로 추정됐다는 사실로부터 $t_j$ 또한 구할 수 있습니다. 이로부터 SGD로 구한 $β$가 얼마나 믿을 만한지 추론해낼 수 있습니다.

정규화

Regulization은 회귀계수에 제약을 가해 일반화 성능을 높이는 기법입니다. 정규화와 관련해서는 이곳을 참고하시면 좋을 것 같습니다. 이 가운데 $β$의 L2 norm을 제한하는 릿지 회귀를 살펴보겠습니다. 파이썬 코드는 다음과 같습니다.

def ridge_penalty(beta, alpha):
    # alpha는 페널티의 강도를 조절하는 하이퍼 파라메터
    # ridge회귀는 상수에 대한 패널티는 주지 않는다
    return alpha * dot(beta[1:], beta[1:])

def squared_error_ridge(x_i, y_i, beta, alpha):
    # beta를 사용할 때 오류와 페널티의 합을 반환
    return error(beta, x_i, y_i) ** 2 + 
			ridge_penalty(beta, alpha)

def ridge_penalty_gradient(beta, alpha):
    # ridge회귀는 상수에 대한 패널티는 주지 않는다
    return [0] + [2 * alpha * beta_j for beta_j in beta[1:]]

def vector_add(v, w):
    """adds two vectors componentwise"""
    return [v_i + w_i for v_i, w_i in zip(v,w)]

def squared_error_ridge_gradient(x_i, y_i, beta, alpha):
    # i번째 오류 제곱값과 패널티 합의 기울기
    return vector_add(squared_error_gradient(x_i, y_i, beta),
                      ridge_penalty_gradient(beta, alpha))

def estimate_beta_ridge(x, y, alpha):
    # 페널티 파라메터가 alpha인 ridge 회귀를 SGD로 학습
    from functools import partial
    beta_initial = [random.random() for _ in x[0]]
    return minimize_stochastic(partial(squared_error_ridge, 
                                       alpha=alpha),
                               partial(squared_error_ridge_gradient, 
                                       alpha=alpha),
                               x, y,
                               beta_initial,
                               0.001)

alpha=0으로 실행하면 기존 선형회귀와 동일한 결과가 나옵니다. alpha를 증가시킬 수록 $β$가 작아집니다.

random.seed(0)
beta_0 = estimate_beta_ridge(input, output, alpha=0.0)
beta_10 = estimate_beta_ridge(input, output, alpha=10.0)

alpha=10으로 실행하면 $β$는 [28.30817306438882, 0.7309844552226781, -0.9146285684215772, -0.01693934899423252]으로 추정되는데요. 0에 가까운 네번째 회귀계수에 해당하는 변수는 ‘박사학위 취득 여부’입니다. 다시 말해 박사학위 취득 여부는 사이트 이용시간에 큰 영향을 주는 변수가 아니라는 걸 알 수 있습니다.

Comment  Read more

로지스틱회귀 파라메터 추정

02 Jul 2017 | regression

이번 글에서는 로지스틱회귀 모델의 계수를 추정하는 방법을 살펴보도록 하겠습니다. 이번 글은 고려대 김성범 교수님 강의와 ‘밑바닥부터 시작하는 데이터과학(조엘 그루스 지음, 인사이트 펴냄)’을 정리하였음을 먼저 밝힙니다. 로지스틱회귀에 대한 기본적인 내용은 이곳을 참고하시면 좋을 것 같습니다. 그럼 시작하겠습니다.

최대우도추정법

로지스틱회귀 모델 계수를 추정하기 위해서는 우선 최대우도추정법(Maximum Likelihood Estimation) 개념을 먼저 짚고 넘어가야 합니다. 가령 데이터가 임의의 파라메터 $θ$에 의존하는 확률분포를 따르고 표본데이터 $v_1, v_2,…,v_n$이 주어졌다고 해보겠습니다.

우리는 데이터가 발생할 확률 $P(v_1,…,v_n$|$θ)$를 구하고 싶은데, $θ$를 모르는 상황입니다. 따라서 우리는 베이즈 규칙을 활용해 다음과 같이 $θ$가 발생할 우도(likelihood)로 바꿔서 생각할 수 있습니다. 다시 말해 조건절과 결과절을 뒤집어 원하는 값을 구하는 것입니다. \(L(\theta |{ v }_{ 1 },{ v }_{ 2 },...,{ v }_{ n })\) 이 경우 가장 적절한 $θ$는 우도를 최대화해주는 값이 될 것입니다. 다시 말해 관측된 데이터가 발생할 경우를 가장 높게 만들어주는 값이라는 의미입니다.

전제 : 베르누이 분포

로지스틱 회귀는 베르누이 시행(Bernoulli trial)을 전제로 하는 모델입니다. 베르누이 시행이란 어떤 실험이 두 가지 결과만을 가지는 실험을 가리킵니다. 베르누이 시행의 결과에 따라 0(실패) 또는 1(성공)의 값을 대응시키는 확률변수(random variable)를 베르누이 확률변수라 합니다. 이 확률변수의 확률분포를 베르누이 분포라고 합니다. 베르누이 확률변수 Y의 분포는 아래 표와 같습니다.

$Y$ 0 1
$P(Y=y_i)$ $1- p$ $p$

위 표를 수식으로 정리하면 아래와 같습니다.

[P(Y={ y }{ i })={ p }^{ { y }{ i } }{ (1-p) }^{ 1-{ y }{ i } }\quad ({ y }{ i }=0,1)]

베르누이 확률변수 Y에 관한 우도함수(likelihood function)는 아래와 같습니다.

[L=\prod { i }^{ }{ { p }^{ { y }{ i } }{ (1-p) }^{ 1-{ y }_{ i } }}]

로지스틱회귀의 우도함수

학습데이터에 관측치가 $i$개가 있고, 정답 범주가 2개(1 또는 0)뿐인 이항로지스틱 모델의 파라메터 $β$가 주어졌다고 가정해 보겠습니다. 그러면 $i$번째 관측치의 종속변수 $y_i$는 $σ(β^Tx_i)$의 확률로 1, $1-σ(β^Tx_i)$의 확률로 0이 됩니다. 여기에서 $x_i$는 $i$번째 관측치의 독립변수, $σ$는 로지스틱 함수를 가리킵니다. 따라서 로지스틱회귀의 우도함수는 다음과 같이 쓸 수 있습니다.

[L=\prod { i }^{ }{ { \sigma ({ \beta }^{ T }\overrightarrow { { x }{ i } } ) }^{ { y }{ i } }{ \left{ 1-\sigma ({ \beta }^{ T }\overrightarrow { { x }{ i } } ) \right} }^{ 1-{ y }_{ i } } }]

로지스틱 회귀의 파라메터 $β$는 앞서 언급한 MLE로 구합니다. 로그는 단조 증가함수이므로 로그 우도함수(log-likelihood function)를 최대로 하는 회귀계수 $β$는 동시에 우도를 최대화하는 $β$이며 그 역도 성립합니다. 로그 우도함수는 아래와 같이 정리할 수 있습니다.

[\ln { L } =\sum { i }^{ }{ { y }{ i }\ln { \left{ \sigma ({ \beta }^{ T }{ \overrightarrow { x_{ i } } }) \right} } } +\sum { i }^{ }{ \left( 1-{ y }{ i } \right) \ln { \left{ 1-\sigma ({ \beta }^{ T }{ \overrightarrow { x_{ i } } }) \right} } }]

다만 위 로그 우도함수는 추정 대상 파라메터인 회귀계수 $β$에 대해 비선형이기 때문에 선형회귀와 같이 명시적인 해가 존재하지 않습니다. 따라서 Stochastic Gradient Descent(SGD) 같은 반복적이고 점진적인 방식으로 해를 구하게 됩니다.

Stochastic Gradient Descent

SGD는 반복문을 돌 때마다 개별 데이터 포인트에 대한 그래디언트를 계산하고 이 그래디언트의 반대 방향으로 파라메터를 업데이트해 함수의 최소값을 구하는 기법입니다. 로지스틱 회귀의 파라메터는 명시적인 해가 존재하지 않기 때문에 이처럼 값을 조금씩 바꿔가면서 최적해를 찾아가는 것입니다.

그런데 우리는 로지스틱 회귀의 로그 우도를 최대화하는 파라메터 $β$를 찾으려는 것이므로 기존 SGD의 반대 방향으로 업데이트를 수행해야 할 것입니다. SGD 관련 메인 파이썬 코드는 다음과 같습니다.

def minimize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
    # SGD 방식으로 gradient descent
    # minimize_batch보다 훨씬 빠르다

    data = zip(x, y)
    # theta_0를 초기 theta로
    theta = theta_0
    # alpha_0를 초기 이동거리(step_size)로
    alpha = alpha_0
    # 시작할 때의 최소값
    min_theta, min_value = None, float("inf")
    iterations_with_no_improvement = 0

    # 만약 100번 넘게 반복하는 동안 value가 더 작아지지 않으면 멈춤
    while iterations_with_no_improvement < 100:
        value = sum( target_fn(x_i, y_i, theta) for x_i, y_i in data )

        # 새로운 최솟값을 찾았다면
        if value < min_value:
            # 이 값을 저장
            min_theta, min_value = theta, value
            # 100번 카운트도 초기화
            iterations_with_no_improvement = 0
            # 기본 이동거리로 돌아감
            alpha = alpha_0

        # 만약 최솟값이 줄어들지 않는다면
        else:
            # 이동거리 축소
            alpha *= 0.9
            # 100번 카운트에 1을 더함
            iterations_with_no_improvement += 1
		
        # 반복문이 돌 때마다 in_random_order를 호출하기 때문에
        # 매 iter마다 그래디언트를 계산하는 순서가 달라짐
        for x_i, y_i in in_random_order(data):
            # 각 데이터 포인트에 대해 그래디언트를 계산
            gradient_i = gradient_fn(x_i, y_i, theta)
            # 기존 theta에서, 학습률(alpha)과 그래디언트를 뺀 것을 업데이트
            theta = vector_subtract(theta, scalar_multiply(alpha, gradient_i))

    return min_theta

# maximize_stochastic = -minimize_stochastic
def maximize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
    return minimize_stochastic(negate(target_fn),
                               negate_all(gradient_fn),
                               x, y, theta_0, alpha_0)
# x를 입력하면 -f(x)를 반환해주는 함수
def negate(f):
    return lambda *args, **kwargs: -f(*args, **kwargs)

# f가 여러 숫자를 반환할 때 모든 숫자를 음수로 변환
def negate_all(f):
    return lambda *args, **kwargs: [-y for y in f(*args, **kwargs)]

그러면 이제는 target_fn과 gradient_fn, theta를 정의해야 합니다. target_fn은 목적함수로 로그 우도를 아웃풋으로 산출하는 함수를 의미합니다. 로지스틱회귀 모델의 로그 우도함수는 다음과 같습니다.

import math

# 로지스틱 함수
def logistic(x):
    return 1.0 / (1 + math.exp(-x))

# 로지스틱 함수의 1차 도함수
def logistic_prime(x):
    return logistic(x) * (1 - logistic(x))

# i번째 관측치에 대한 로그우도
def logistic_log_likelihood_i(x_i, y_i, beta):
    # y_i의 레이블이 1이면
    if y_i == 1:
        return math.log(logistic(dot(x_i, beta)))
    # y_i의 레이블이 0이면
    else:
        return math.log(1 - logistic(dot(x_i, beta)))

# 데이터 전체에 대한 로그우도 (전체 합)
def logistic_log_likelihood(x, y, beta):
    return sum(logistic_log_likelihood_i(x_i, y_i, beta)
               for x_i, y_i in zip(x, y))

‘gradient_fn’은 각 파라메터에 대한 목적함수의 그래디언트를 가리킵니다. 미분 공식을 활용해 로그우도 함수를 $β$에 대해 해석적으로 미분한 식을 코드로 표현한 것입니다.

def logistic_log_partial_ij(x_i, y_i, beta, j):
    """here i is the index of the data point,
    j the index of the derivative"""

    return (y_i - logistic(dot(x_i, beta))) * x_i[j]
    
def logistic_log_gradient_i(x_i, y_i, beta):
    """the gradient of the log likelihood 
    corresponding to the i-th data point"""

    return [logistic_log_partial_ij(x_i, y_i, beta, j)
            for j, _ in enumerate(beta)]
            
def logistic_log_gradient(x, y, beta):
    return reduce(vector_add,
                  [logistic_log_gradient_i(x_i, y_i, beta)
                   for x_i, y_i in zip(x,y)])

마지막으로 ‘theta’는 해당 목적함수의 파라메터인데요, 우리 문제에선 $β$를 말합니다. 만약 독립변수가 2개라면 상수항을 포함해 3차원 벡터가 될 것입니다. 이밖에 위 코드 실행에 필요한 함수는 다음과 같습니다.

def in_random_order(data):
    # Stochastic Gradient Descent 수행을 위한 함수로,
    # 한번 반복문을 돌 때마다 임의의 순서로 데이터 포인트를 반환
    # 데이터 포인트의 인덱스를 list로 생성
    indexes = [i for i, _ in enumerate(data)]
    # 이 인덱스를 랜덤하게 섞는다
    random.shuffle(indexes)
    # 이 순서대로 데이터를 반환한다
    for i in indexes:
        yield data[i]

def vector_subtract(v, w):
    """subtracts two vectors componentwise"""
    return [v_i - w_i for v_i, w_i in zip(v,w)]

def scalar_multiply(c, v):
    return [c * v_i for v_i in v]

def dot(v, w):
    """v_1 * w_1 + ... + v_n * w_n"""
    return sum(v_i * w_i for v_i, w_i in zip(v, w))

분석 예시

분석 대상 데이터는 다음과 같습니다. 독립변수는 2개, 독립변수는 범주형 1개(0 또는 1)입니다.

data = [(0.7,48000,1),(1.9,48000,0),(2.5,60000,1),(4.2,63000,0),(6,76000,0),(6.5,69000,0),(7.5,76000,0),(8.1,88000,0),(8.7,83000,1),(10,83000,1),(0.8,43000,0),(1.8,60000,0),(10,79000,1),(6.1,76000,0),(1.4,50000,0),(9.1,92000,0),(5.8,75000,0),(5.2,69000,0),(1,56000,0),(6,67000,0),(4.9,74000,0),(6.4,63000,1),(6.2,82000,0),(3.3,58000,0),(9.3,90000,1),(5.5,57000,1),(9.1,102000,0),(2.4,54000,0),(8.2,65000,1),(5.3,82000,0),(9.8,107000,0),(1.8,64000,0),(0.6,46000,1),(0.8,48000,0),(8.6,84000,1),(0.6,45000,0),(0.5,30000,1),(7.3,89000,0),(2.5,48000,1),(5.6,76000,0),(7.4,77000,0),(2.7,56000,0),(0.7,48000,0),(1.2,42000,0),(0.2,32000,1),(4.7,56000,1),(2.8,44000,1),(7.6,78000,0),(1.1,63000,0),(8,79000,1),(2.7,56000,0),(6,52000,1),(4.6,56000,0),(2.5,51000,0),(5.7,71000,0),(2.9,65000,0),(1.1,33000,1),(3,62000,0),(4,71000,0),(2.4,61000,0),(7.5,75000,0),(9.7,81000,1),(3.2,62000,0),(7.9,88000,0),(4.7,44000,1),(2.5,55000,0),(1.6,41000,0),(6.7,64000,1),(6.9,66000,1),(7.9,78000,1),(8.1,102000,0),(5.3,48000,1),(8.5,66000,1),(0.2,56000,0),(6,69000,0),(7.5,77000,0),(8,86000,0),(4.4,68000,0),(4.9,75000,0),(1.5,60000,0),(2.2,50000,0),(3.4,49000,1),(4.2,70000,0),(7.7,98000,0),(8.2,85000,0),(5.4,88000,0),(0.1,46000,0),(1.5,37000,0),(6.3,86000,0),(3.7,57000,0),(8.4,85000,0),(2,42000,0),(5.8,69000,1),(2.7,64000,0),(3.1,63000,0),(1.9,48000,0),(10,72000,1),(0.2,45000,0),(8.6,95000,0),(1.5,64000,0),(9.8,95000,0),(5.3,65000,0),(7.5,80000,0),(9.9,91000,0),(9.7,50000,1),(2.8,68000,0),(3.6,58000,0),(3.9,74000,0),(4.4,76000,0),(2.5,49000,0),(7.2,81000,0),(5.2,60000,1),(2.4,62000,0),(8.9,94000,0),(2.4,63000,0),(6.8,69000,1),(6.5,77000,0),(7,86000,0),(9.4,94000,0),(7.8,72000,1),(0.2,53000,0),(10,97000,0),(5.5,65000,0),(7.7,71000,1),(8.1,66000,1),(9.8,91000,0),(8,84000,0),(2.7,55000,0),(2.8,62000,0),(9.4,79000,0),(2.5,57000,0),(7.4,70000,1),(2.1,47000,0),(5.3,62000,1),(6.3,79000,0),(6.8,58000,1),(5.7,80000,0),(2.2,61000,0),(4.8,62000,0),(3.7,64000,0),(4.1,85000,0),(2.3,51000,0),(3.5,58000,0),(0.9,43000,0),(0.9,54000,0),(4.5,74000,0),(6.5,55000,1),(4.1,41000,1),(7.1,73000,0),(1.1,66000,0),(9.1,81000,1),(8,69000,1),(7.3,72000,1),(3.3,50000,0),(3.9,58000,0),(2.6,49000,0),(1.6,78000,0),(0.7,56000,0),(2.1,36000,1),(7.5,90000,0),(4.8,59000,1),(8.9,95000,0),(6.2,72000,0),(6.3,63000,0),(9.1,100000,0),(7.3,61000,1),(5.6,74000,0),(0.5,66000,0),(1.1,59000,0),(5.1,61000,0),(6.2,70000,0),(6.6,56000,1),(6.3,76000,0),(6.5,78000,0),(5.1,59000,0),(9.5,74000,1),(4.5,64000,0),(2,54000,0),(1,52000,0),(4,69000,0),(6.5,76000,0),(3,60000,0),(4.5,63000,0),(7.8,70000,0),(3.9,60000,1),(0.8,51000,0),(4.2,78000,0),(1.1,54000,0),(6.2,60000,0),(2.9,59000,0),(2.1,52000,0),(8.2,87000,0),(4.8,73000,0),(2.2,42000,1),(9.1,98000,0),(6.5,84000,0),(6.9,73000,0),(5.1,72000,0),(9.1,69000,1),(9.8,79000,1),]

전체 코드를 구동하는 명령은 다음과 같습니다.

import random
data = map(list, data)  # change tuples to lists
x = [[1] + row[:2] for row in data]  # each element is [1, experience, salary]
y = [row[2] for row in data]  # each element is paid_account
beta_0 = [1, 1, 1]
beta_hat = maximize_stochastic(logistic_log_likelihood_i,
                                logistic_log_gradient_i,
                                x, y, beta_0)

Comment  Read more

베이즈 규칙과 다양한 문제들

01 Jul 2017 | Bayes

이번 글에서는 베이즈 규칙을 활용해 풀 수 있는 여러 문제를 살펴보도록 하겠습니다. 이 글은 고려대 김성범 교수님 강의와 ‘Think Bayes(앨런 B. 다우니 지음, 권정민 옮김, 한빛미디어 펴냄)’를 정리했음을 먼저 밝힙니다. 베이지안 추론과 관련해서는 이곳과 이 글, 문서 분류를 위한 나이브 베이지안 분류기와 관련해서는 이곳을 참고하시면 좋을 것 같습니다. 그럼 시작하겠습니다.

베이즈 규칙

일반적으로 사건 $A_1$, $A_2$, $A_3$가 서로 배반(mutually exclusive)이고 $A_1$, $A_2$, $A_3$의 합집합이 표본공간(sample space)과 같으면 사건 $A_1$, $A_2$, $A_3$는 표본공간 $S$의 분할이라고 정의합니다. 우리가 관심있는 사건 $B$가 나타날 확률을 그림과 식으로 나타내면 다음과 같습니다.

[P(B)=P({ A }{ 1 }\cap B)+P({ A }{ 2 }\cap B)+P({ A }_{ 3 }\cap B)]

$P(B)$를 조건부확률의 정의를 이용해 다시 쓰면 아래와 같습니다. 이를 전확률 공식(Law of Total Probability) 또는 베이즈 법칙이라고 합니다.

[P(B)=P({ A }_{ 1 })P(B { A }{ 1 })+P({ A }{ 2 })P(B { A }{ 2 })+P({ A }{ 3 })P(B { A }{ 3 })=\sum _{ i=1 }^{ 3 }{ P({ A }{ i })P(B { A }_{ i }) }]

보통 $P(A_1)$, $P(A_2)$, $P(A_3)$는 미리 알고 있다는 의미의 사전확률(prior probability)로 불립니다. $P(B$|$A1)$, $P(B$|$A2)$, $P(B$|$A3)$는 우도(likelihood probability)라 부릅니다.

그럼 우리가 관심있는 사건인 $B$가 $A_1$에 기인했을 조건부확률은 어떻게 구할까요? 바로 아래와 같이 구할 수 있습니다.

[\begin{align} P({ A }_{ 1 }|B)&=\frac { P({ A }_{ 1 })P(B|{ A }_{ 1 }) }{ P(B) }
&=\frac { P({ A }_{ 1 })P(B|{ A }_{ 1 }) }{ P({ A }_{ 1 })P(B|{ A }_{ 1 })+P({ A }_{ 2 })P(B|{ A }_{ 2 })+P({ A }_{ 3 })P(B|{ A }_{ 3 }) } \end{align}]

$P(A1$|$B)$는 사건 B를 관측한 후에 그 원인이 되는 사건 $A$의 확률을 따졌다는 의미의 사후확률(posterior probability)로 정의됩니다. 사후확률은 사건 $B$의 정보가 더해진, 사전확률의 업데이트 버전 정도라고 생각하면 좋을 것 같습니다. (Posterior probability is an updated version of prior probability) 같은 방식으로 $P(A2$|$B)$, $P(A3$|$B)$도 구할 수 있습니다.

한번 예를 들어보겠습니다. 어떤 이름모를 질병에 걸린 환자가 전체 인구의 약 1% 정도 되는 것으로 알려져 있다고 칩시다. 그렇다면 전체 인구라는 표본공간에서 질병에 걸릴 확률 P(D)는 0.01, 그렇지 않을 확률 P(~D)는 0.99입니다. 이것이 바로 우리가 이미 알고 있는 사전확률이 되겠네요.

질병 발생 여부를 측정해주는 테스트의 정확도는 이렇다고 합니다. 진짜 환자를 양성이라고 정확하게 진단할 확률 P(+|D)은 97%, 정상환자를 양성이라고 오진할 확률 P(+|~D)은 6%입니다. 이것이 바로 우도입니다.

그럼 진단테스트 결과 양성이라고 나왔는데 실제 환자일 확률 P(D|+)는 얼마일까요? 이건 우리가 이미 알고 있는 정보를 활용해 아래와 같이 구할 수 있습니다. 이것이 바로 사후확률입니다. 즉 ‘진단’이라는 사건의 정보를 더해 ‘질병에 걸릴 확률’이라는 사전확률을 업데이트한거죠.

[\begin{align} P(+)&=P(D\cap +)+P(\sim D\cap +)\ &=P(D)P(+|D)+P(\sim D)P(+|\sim D)\ &=0.01\times 0.97+0.99\times 0.06\ &=0.691\ P(D|+)&=\frac { P(D)P(+|D) }{ P(+) }
&=\frac { 0.01\times 0.97 }{ 0.691 }
&=0.014 \end{align}]

그런데 왜 이렇게 복잡하게 사후확률을 구하는거냐고요? 실제로 사후확률은 구하기 어려운 경우가 많다고 합니다. 그에 반해 사전확률은 우리가 이미 알고 있는 값이고, 우도는 비교적 계산하기 수월합니다. 그래서 전확률법칙과 사전확률과 우도를 활용해 사후확률을 도출해내는 것이죠.

다시 말해 우리가 알고 싶은 확률을 단박에 계산하기가 까다로울 때 조건과 결과를 뒤집어서 우회적으로 계산하는 것, 이것이 베이즈 모델의 강점이라고 말할 수 있겠습니다.

쿠키 문제

쿠키가 들어 있는 그릇 두 개가 있다고 가정해보겠습니다. 첫번째 그릇에는 바닐라 쿠키 30개와 초콜렛 쿠키 10개가 들어있고, 두번째 그릇에는 두 가지 쿠키가 종류별로 20개씩 들어 있습니다.

어떤 그릇인지 보지 않고 한 그릇에서 임의로 쿠키를 집었는데 바닐라 쿠키였다고 칩시다. 그렇다면 이 때 ‘이 바닐라 쿠키가 그릇 1에서 나왔을 가능성’은 얼마일까요? 사실 이 확률을 계산하기는 쉽지 않습니다.

그런데 ‘그릇1에서 바닐라 쿠키가 나올 확률’은 30/40으로 구하기 쉽습니다. 베이즈 규칙을 활용해 조건절과 결과절을 바꾸어 원하는 확률을 구해보겠습니다. 저 아래 표의 첫번째 열을 기준으로 계산한 식과 용어 설명은 다음과 같습니다.

[P(H D)=\frac { P(H)P(D H) }{ P(D) } =\frac { \frac { 1 }{ 2 } \times \frac { 3 }{ 4 } }{ \frac { 5 }{ 8 } } =\frac { 3 }{ 5 }]

$P(H)$ : 어떤 쿠키를 골랐던지 상관없이 그릇1을 골랐을 확률. 문제에서는 그릇을 임의로 선택한 것이므로 0.5라고 가정할 수 있습니다. 이를 데이터를 보기 전의 가설의 확률, 즉 사전확률입니다.

$P(D$|$H)$ : 그릇1에서 바닐라 쿠키가 나올 확률. 3/4입니다. 이를 데이터가 가설에 포함될 확률, 즉 우도입니다.

$P(D)$ : 바닐라 쿠키를 고를 확률입니다. 각 그릇을 고를 확률이 동일하고 그릇에 동일한 쿠키가 들어있으므로 어떤 그릇을 택하든 바닐라 쿠키를 고를 확률도 같습니다. 두 그릇에 50개 바닐라 쿠키와 30개의 초콜렛 쿠키가 들어있으므로 $P(V)$는 5/8이 됩니다. 이를 어떤 가설에든 포함되는 데이터의 비율, 즉 한정상수입니다.

$P(H$|$D)$ : 바닐라 쿠키가 그릇1에서 나왔을 확률. 우리가 알고 싶은 확률입니다. 이를 데이터를 확인한 이후의 가설 확률, 즉 사후확률입니다. 이는 베이즈 이론을 통시적(diachronic)으로 해석한 것으로 가설에 대한 확률이 시간에 따라 새로운 데이터를 접하게 되면서 달라진다, 다시 말해 데이터 $D$의 관점에서 봤을 때 가설 $H$의 확률을 업데이트해 준다는 방식으로 이해하는 것입니다.

여기에서 한정상수 $P(D)$는 다음과 같은 두 가지 조건을 만족할 경우 어떤 가설이든지 같은 값을 지니게 돼 계산을 생략할 수 있습니다.

상호배제(mutually exclusive) : 집합 중 하나의 가설만 참이다

전체포괄(collectively exhaustive) : 고려 대상 가설 이외에 다른 가설이 전혀 없는 경우

쿠키 문제를 가설별로 표로 나타내면 다음과 같습니다.

항목 가설1(그릇1) 가설2(그릇2)
사전확률 $P(H)$ 1/2 1/2
우도 $P(D$|$H)$ 3/4 1/2
사전확률 × 우도 3/8 1/4
한정상수 $P(D)$ 5/8 5/8
사후확률 $P(H$|$D)$ 3/5 2/5

위 표에서 한정상수 $P(D)$ 계산을 생략할 수 있는 이유에 대해 살펴보겠습니다. 사후확률은 사전확률에 우도를 곱한 뒤 한정상수를 나누어 구합니다. 그런데 한정상수가 가설에 상관없이 같은 값을 가진다면 사전확률에 우도를 곱한 걸 모두 더한 값이 1이 되도록 정규화하면 한정상수 계산을 생략해도 동일한 결과가 나옵니다.

위 표에서 사전확률에 우도를 곱한 값은 그릇1이 3/8, 그릇2가 1/4입니다. 이를 모두 더하면 5/8이 되는데요. 이 값이 1이 되도록 정규화를 하려면 그릇1과 그릇2가 가진 값에 8/5를 각각 곱해주면 됩니다. 그런데 이 값은 정확히 사후확률과 일치합니다. 한정상수가 가설에 관계없이 동일하다면 계산을 생략해도 관계 없다는 이야기입니다.

쿠키문제를 파이썬 코드로 풀어보겠습니다.

import thinkbayes as tb

class Cookie(tb.Pmf):
    """A map from string bowl ID to probablity."""
	
    # 사전확률 정의
    def __init__(self, hypos):
        """Initialize self.

        hypos: sequence of string bowl IDs
        """
        tb.Pmf.__init__(self)
        for hypo in hypos:
            # 가설마다 동일한 사전확률(1) 부여
            self.Set(hypo, 1)
        # 총합이 1이 되도록(확률이 되도록) 정규화
        # 가설이 두 개이므로 정규화 결과 가설마다 0.5 사전확률 가짐
        self.Normalize()

    # 사후확률 구하기
    def Update(self, data):
        """Updates the PMF with new data.

        data: string cookie type
        """
        for hypo in self.Values():
            # 우도 구하기
            like = self.Likelihood(data, hypo)
            # 사전확률에 우도 업데이트(곱)
            self.Mult(hypo, like)
        # 사전확률 x 우도 모든 값의 합이 1이 되도록 정규화
        # 한정상수가 동일하므로 사후확률 계산에서 생략
        self.Normalize()

    mixes = {
        'Bowl 1':dict(vanilla=0.75, chocolate=0.25),
        'Bowl 2':dict(vanilla=0.5, chocolate=0.5),
        }

    # 우도 함수
    def Likelihood(self, data, hypo):
        """The likelihood of the data under the hypothesis.

        data: string cookie type
        hypo: string bowl ID
        """
        # mixes 데이터 자체가 우도를 나타내고 있으므로 
        # 그냥 불러오기만 하면 우도함수 정의 끝
        mix = self.mixes[hypo]
        like = mix[data]
        return like

# 가설정의 (가설1=그릇1, 가설2=그릇2)
hypos = ['Bowl 1', 'Bowl 2']

# 가설에 따른 객체 선언
pmf = Cookie(hypos)

# 사후확률 구하기
# 사전확률에 우도 업데이트 + 정규화
pmf.Update('vanilla')

# 결과 출력
for hypo, prob in pmf.Items():
    print hypo, prob

M&M 문제

여러 색의 설탕이 입혀진 작은 초콜렛 M&M은 만드는 회사는 시간에 따라 색 조합을 바꿔 왔습니다. 이 회사는 1995년에 파란색을 추가했는데요. 그 전에는 봉지 내 색 조합이 갈색 30%, 노랑 20%, 빨강 20%, 녹색 10%, 주황 10%, 황갈색 10%였습니다. 파란색을 추가한 뒤에는 파랑 24%, 녹색 20%, 주황 16%, 노랑 14%, 빨강 13%, 갈색 13%가 됐습니다.

한 친구가 M&M을 두 봉지 샀는데 각각 생산년도가 1994년, 1996년이라고 칩시다. 생산년도를 알려주지 않고 각 봉지에서 M&M을 하나씩 꺼냈을 때 한 알은 노란색이고 한 알은 녹색이었습니다. 문제입니다. 이 때 노랑 초콜렛이 1994년에 생산한 봉지에서 나왔을 확률은 얼마일까요?

먼저 가설을 만들어보겠습니다. 노란 초콜렛은 봉지1에서 녹색 초콜렛은 봉지2에서 꺼냈다고 칩시다. 이 때 가설은 다음과 같습니다.

가설1 : 봉지1은 1994년에 생산했고 봉지2는 1996년에 생산했다.

가설2 : 봉지1은 1996년에 생산했고 봉지2는 1994년에 생산했다.

이를 쿠키문제와 마찬가지로 표로 만들면 다음과 같습니다. 쿠키문제와 마찬가지로 한정상수가 같은 값이기 때문에 계산을 생략할 수 있습니다.

항목 가설1 가설2
사전확률 $P(H)$ 1/2 1/2
우도 $P(D$|$H)$ 0.2×0.2 0.14×0.1
사전확률 × 우도 0.02 0.007
사후확률 $P(H$|$D)$ 20/27 7/27

M&M 문제의 파이썬 풀이는 다음과 같습니다. Think Bayes 저자는 베이즈 규칙을 언제든 사용할 수 있도록 Suite 객체로 캡슐화하였습니다. 쿠키문제의 풀이와 본질적으로 다르지 않다는 얘기입니다. 어쨌든 Suite 객체로 베이즈 문제를 풀 때는 데이터의 우도 값과 우도 함수만 정의해 주면 됩니다.

from thinkbayes import Suite


class M_and_M(Suite):
    """Map from hypothesis (A or B) 
    to probability."""
	
    # 1994년 데이터
    mix94 = dict(brown=30,
                 yellow=20,
                 red=20,
                 green=10,
                 orange=10,
                 tan=10)
	
    # 1996년 데이터
    mix96 = dict(blue=24,
                 green=20,
                 orange=16,
                 yellow=14,
                 red=13,
                 brown=13)
	
    # 가설 정의
    hypoA = dict(bag1=mix94, bag2=mix96)
    hypoB = dict(bag1=mix96, bag2=mix94)
    hypotheses = dict(A=hypoA, B=hypoB)
	
    # 우도함수 정의
    def Likelihood(self, data, hypo):
        """Computes the likelihood of the data 
        under the hypothesis.
        hypo: string hypothesis (A or B)
        data: tuple of string bag, string color
        """
        bag, color = data
        # 데이터가 우도이므로 참조하기만 하면 됨
        mix = self.hypotheses[hypo][bag]
        like = mix[color]
        return like

# 가설과 객체 선언
suite = M_and_M('AB')

# 사전확률에 데이터를 반영해 사후확률 계산
suite.Update(('bag1', 'yellow'))
suite.Update(('bag2', 'green'))

# 결과 출력
suite.Print()

주사위 문제

4면체, 6면체, 8면체, 12면체, 20면체 주사위가 든 상자가 있다고 가정해 봅시다. 상자에서 임의로 주사위 하나를 집어서 던졌더니 12가 나왔습니다. 그러면 각 주사위를 선택했을 확률은 어떻게 될까요?

여기에서 사전확률은 주사위를 선택할 확률, 우도는 해당 주사위를 던졌을 때 12가 나올 확률입니다. 이를 표로 나타내면 다음과 같습니다.

항목 4면체 6면체 8면체 12면체 20면체
사전확률 $P(H)$ 1/5 1/5 1/5 1/5 1/5
우도 $P(D$|$H)$ 0 0 0 1/12 1/20
사전확률 × 우도 0 0 0 1/60 1/100
사후확률 $P(H$|$D)$ 0 0 0 5/8 3/8

파이썬 코드는 다음과 같습니다.

import thinkbayes as tb

class Dice(tb.Suite):
    # 우도함수 정의
    def Likelihood(self, data, hypo):
        if hypo < data:
            return 0
        else:
            return 1.0/hypo

# 가설정의 및 객체 선언
suite=Dice([4,6,8,12,20])

# 사전확률에 우도 곱해줌
# 우도는 관찰결과(데이터)
# 이를 통해 사후확률 계산
suite.Update(12)

# 결과 출력
suite.Print()

위 코드의 실행 결과는 다음과 같습니다.

4 0.0 6 0.0 8 0.0 12 0.625 20 0.375

기본코드

객체 생성에 쓰인 코드는 thinkbayes.py인데 정리 겸 하단에 첨부하였습니다. 출처는 이곳입니다.

Comment  Read more