HAMA 블로그

데이터 분석에 대한 기본적인 감과 덤으로 파이썬 코딩에 대한 감을 익히기 위한 강좌입니다. 
개발자는 코드로 이해하는게 가장 빠른 길이며, 오래 기억하는 길이라 믿습니다.

DTW (Dynamic time wrapping) 란?

두개의 시계열 데이터가 있다고 할 때 그 둘간의 유사도를 알아내기 위한 알고리즘 중 하나 이다.

예를들어 IoT 센서에서 시간마다 생성되는 데이터를 비교하여, 현재 데이터흐름은 과거의 어떤 흐름과 매칭되는지 체크하는데 사용된다. 즉 에어컨을 켰을 때 전력패턴을 기억해서, 후에 내 스마트폰을 통해 외부에서 전력사용을 감지해서 누가 내 집에서 에어컨을 켠다던가 하는?? 도둑놈?? 혹은 냉장고 문이 열렸을 때의 패턴을 감지하여  혹시 아들내미가 냉장고 문을 열어두고 가진 않았는지? (물론 DTW 만 가지고서 저런 NILM 문제를 해결할 수 있는것은 아니다. ) 가장 많이 떠올릴 수 있는 데이터인 주식데이터에도 사용 될 수 있다.

다음 그림을 보자. 

보통 떠올리기에는 0~1분 사이의 1초간의 데이터가 있다고 할 때, 0초,1초,2초 각각의 데이터의 차의 제곱의 평균  를 이용해서 차이를 구하는게 가장 먼저 떠오를 것이다. 하지만 그렇게 되면 위 그림의 오른쪽과 같은 좀 비틀어진 상황에서는 영 맥을 못추는 패턴매칭이 될 것이다.DTW 는 이럴때 사용되는 시계열 패턴 매칭 알고리즘이다. 어떻게 해결 할 수 있을지 잠시 생각을 해보자.

자신과 같은 시간 index 를 가진 요소와의 비교 뿐 만 아니라, 그 주변의 다른 요소와도 비교를 해서 더 비슷한 요소를 자신의 비교쌍으로 놓는다면 어느정도 극복되지 않을까? 자세한것은 코드를 보고 이해해보자.

사실 쉬운 알고리즘이라도 글로 보면 이해하기가 쉽지 않다. 

코드로 말해요  

# -*- coding: utf-8 -*-

from math import *
import numpy as np
import sys


def DTW(A, B, window=sys.maxint, d=lambda x, y: abs(x - y)):
    # 비용 행렬 초기화 
    A, B = np.array(A), np.array(B)
    M, N = len(A), len(B)
    cost = sys.maxint * np.ones((M, N))

    # 첫번째 로우,컬럼 채우기
    cost[0, 0] = d(A[0], B[0])
    for i in range(1, M):
        cost[i, 0] = cost[i - 1, 0] + d(A[i], B[0])

    for j in range(1, N):
        cost[0, j] = cost[0, j - 1] + d(A[0], B[j])
    # 나머지 행렬 채우기 
    for i in range(1, M):
        for j in range(max(1, i - window), min(N, i + window)):
            choices = cost[i - 1, j - 1], cost[i, j - 1], cost[i - 1, j]
            cost[i, j] = min(choices) + d(A[i], B[j])

    # 최적 경로 구하기 
    n, m = N - 1, M - 1
    path = []

    while (m, n) != (0, 0):
        path.append((m, n))
        m, n = min((m - 1, n), (m, n - 1), (m - 1, n - 1), key=lambda x: cost[x[0], x[1]])

    path.append((0, 0))
    return cost[-1, -1], path


def main():
    A = np.array([1,2,3,4,2,3])
    B = np.array([7,8,5,9,11,9])
 
    cost, path = DTW(A, B, window = 6)
    print 'Total Distance is ', cost
    import matplotlib.pyplot as plt
    offset = 5
    plt.xlim([-1, max(len(A), len(B)) + 1])
    plt.plot(A)
    plt.plot(B + offset)
    for (x1, x2) in path:
        plt.plot([x1, x2], [A[x1], B[x2] + offset])
    plt.show()

if __name__ == '__main__':
    main()

소스참고: https://gist.github.com/bistaumanga/6023705

설명

아래와 같은 2개의 타임시리즈가 있다고 하자.

1. 비용 행렬 초기화 (6x6의 행렬을 만들고, max_int 값으로 채운다)

2 첫번째 로우와 컬럼을 각 타임스탬프(A,B)의 값을 이용하여 채운다.

 -> (0,0) 행렬은 A , B 요소의 첫번째 값의 차로 
 -> A 의 요소와 B의 첫번째 요소간의 차이의 절대값에 이전 값을 더해서 로우를, 그 반대로 컬럼을

3.행렬의 나머지 부분을 채운다.

 -> choices : cost[i - 1, j - 1], cost[i, j - 1], cost[i - 1, j]  이 3요소를 선택한다. (자신을 둘러싼 3방향요소)
 -> choices 에서 선택된 3가지 요소중 작은것과  A[i], B[j] 의 차의 절대값을 더해준다. 

여기서 34라는 DTW 라는 값이 구해졌다.
위에서 B를 np.array([1,3,3,4,3,3]) 이런것으로 한다면 DTW 가 2밖에 되지 않을 것이다.

4. 가장 최적의 경로를 구한다.

 -> 마지막부터 시작해서 주위의 가장 작은 숫자의 칸으로 이동한다.

위의 경우처럼 대각선으로 이동한다는 것은 A,B 두 시계열 데이터가 일정하게 변화했다는 의미이다.

경로를 차트로 그리면 다음과 같다.

B를 np.array([1,3,3,4,3,3]) 이런것으로 한다면 차트는 다음과 같을 것이다.

데이터 분석에 대한 기본적인 감과 덤으로 파이썬 코딩에 대한 감을 익히기 위한 강좌입니다. 
개발자는 코드로 이해하는게 가장 빠른 길이며, 오래 기억하는 길이라 믿습니다.
이 글의 소스인  "밑바닥부터 배우는 데이터과학" 및 "밑바닥부터 시작하는 딥러닝" 은  그런 기본적인 코드를 주요 매개체로 머신러닝을 설명해주는 훌륭한 서적이니 데이터사이언스에 관심있는 개발자라면 추천해드립니다. (수학도 미분,편미분 정도만 알면 되며, 바닥부터 쉽게 설명해줌) 

 * 참고로 "밑바닥부터 배우는 데이터과학" 서적은 numpy,scikit-learn 등의 외부라이브러리를 활용은 배제되었습니다.

k-NN 이란?

k-NN 알고리즘은 가장 쉬운 머신러닝 알고리즘 중의 하나이다. "분류"를 하기위한 알고리즘인데 

예를들어보면 

(1,축구), (2,축구), (3,축구) (7,야구),(8,야구)(10,배구) 이런 데이터 목록들이 있다고 할때 

내가 가지고 있는 데이터인 "3번" 은 축구,야구,배구 중 어디로 분류 될까?  

당연히 축구로 분류 될 것이다.

그럼 4번은 어디로 분류 될까? 4는 3과 가깝고 3번이 축구니깐 또 축구 일 것이다.

축하한다. 당신은 이제 k-NN 알고리즘을 이해했다.!!

그렇다 k-NN 알고리즘은 가까운 거리에 가장 많이 분포하는 것으로 나의 데이터를 분류 하는 알고리즘이다. 따라서 k-NN 알고리즘은 

- 거리를 재는 방법
- 서로 가까운 점들은 유사하다는 가정
- 어느 정도 가까운 거리로 할 것인가? 

정도만 필요하다.

이미지를 통해서 다시 확실히 이해해보자.

자 분류하고 싶은 데이터는 연두색으로 아직 미정이다.
이 데이터는 파랑사각형과 빨강삼각형 중 어디로 분류될까?
자신과 가장 가까운 거리에 있는 것으로 분류를 하면 파랑 사각형으로 분류 될 것이고 
거리를 k=3 으로 정한다면 k=3 이내에 가장 많은 있는 빨강 삼각형으로 분류 될 것이다.
만약  k=3 거리에 동일한 숫자의 사각형,삼각형이 있다면, 거리를 k=2.5 로 접혀서 다시 계산하면 둘 중 하나가 더 많아 질것이다. 혹시 그래도 같다면 k=2 로 줄이고~ 그러는 것이쥐~

 코드로 말해요.

함수1)

def raw_majority_vote(labels):
    votes = Counter(labels)
    winner,_ = votes.most_common(1)[0]
    return winner

위 코드는 데이터목록들을 받아서 , 가장 많이 존재하는 이름을 리턴해준다.
즉 (1,축구),(2,축구)(3,야구) 이런 데이터목록을 받으면 축구를 리턴해준다.

함수2)

def majority_vote(labels):
    vote_counts = Counter(labels):
    winner, winner_count = vote_counts.most_common(1)[0]
    num_winners = len([count 
                                        for count in vote_counts.values() 
                                        if count == winner_count])
    if num_winners == 1:
        return winner
    else:
        return majority_vote(labels[:-1])

위 코드는 데이터목록들을 받아서 , 가장 많이 존재하는 이름을 찾는데, 하나 뿐이라면 그것을 리턴해주고, 
만약 축구,야구가 둘 다 동일한 숫자가 k거리 이내에 있다면 , 제일 멀리 있는 놈을 제외하고 다시 찾는다.

함수3)

* 마지막으로 위의 코드를 이용하여 k-NN 알고리즘을 만들자.

def knn_classify(k, lebeled_points, your_point):
    by_distance = sorted(lebeled_points, 
                                        key=lambda (point, _): distance(point, your_point()))
    k_nearest_labels = [label for _, label in by_distance[:k]]
    
    return majority_vote(k_nearest_labels)

매개변수) 

k : 어느정도 가까운 것들을 찾는가?
labeled_points : 분류에 사용 될 데이터목록들 
your_point : 분류하고 싶은 데이터

설명)

1. 분류에 사용될 데이터들을 분류 될 데이터와 거리 순으로 정렬한다.
2. 정렬된 데이터 중에서 k 거리 이내에 있는 데이터 목록만 따로 majority_vote 에 넘겨서 k 거리이내의 데이터들 중에 가장 많이 포함되 있는 라벨(야구 혹은 축구?) 을 찾는다.

이렇게 당신은 k-NN 알고리즘을 이해하였으며, 구현까지 할 수 있게 되었다.
자신감이 붙었다면 다른 "머신러닝" 알고리즘의 세계로 진출해보자~
다음으로 쉬운 알고리즘으로는 K-Means 를 추천해드린다.

HMM  (은닉마코프법) 

이번 글은  HMM에 대해서 대략 알고 있는 상태에서 , HMM 의 3가지 문제중 가장 고난이도인 학습문제에 대해 코드로 알아보고자 한다. 

HMM 자체에 대해서 알고 싶은 분은 , 오일석 저자의 "패턴인식" 을 추천한다. 이 책은 진짜 설명을 해주려고 애쓰는 모습이 느껴지는 책이고,그 책으로 어느정도 이해한 후에 그에 파생되는 연구(MHMM, FHMM, CFHSMM 등) 에 대한 논문들을 살펴보면 될거 같다.

아래 링크의 글 또한 매우 자세히 HMM 에 대해서 설명하는 글이니 참고하시라~

http://article.sciencepublishinggroup.com/html/10.11648.j.acm.s.2017060401.12.html#paper-content-4 

HMM learning  by code

N (상태후보개수) 과 M (관측후보개수) 과  T 개의 관측열은 고정 되있다고 가정한다.

O = (O0, O1, . . . , OT −1)

1. 초기화

행렬 A,B,π  에 대해서 각각 초기화 한다. π 는 1 × N 의 갯수를 가지고,  A = {aij}

는  N × N 그리고  B = {bj (k)} 는 N × M 을 가진다. 모든 행렬은 row-stochastic 이다.  

참고로 

π  : 초기상태확률에 대한 행렬 (어떤 상태에서 시작될 것인가?) 
A  : 상태전이확률에 대한 행렬 (a 상태에서 b 상태로 이동하는 확률은?)
B  : 관측확률에 대한 행렬 ( a 상태에서 o1,o2,o3 각각이 관측될 확률은?) 

πi ≈ 1/N   즉 초기값 하나하나들은 대략 1/N 의 확률을 갖게 된다. 

aij ≈ 1/N   즉 상태전이확률을 대략 1/N 의 확률을 갖게 된다. 상태가 N 개니깐, N 개의 상태로 이동한다는 말.

bj (k) ≈ 1/M. 즉 관측확률은 대략 1/M 이된다.  M 은 관측후보개수

2. 슈도코드

2. The α-pass

#  α0(i) : 첫번째 관측 (모든 상태)에 대한 전방계산.
#  c0 : 모든 상태에 대한 첫번째 전방변수(확률들의 합) 
c0 = 0
for i = 0 to N − 1   # 모든 상태에 대해서 
    α0(i) = πibi(O0) # 초기상태확률 * 초기상태에서 첫번째 관측이 나올 확률
    c0 = c0 + α0(i)  # 모든 α0(0),α0(1) 의 합  
next i

# α0(i) 에 대한 스케일 조정 
c0 = 1/c0
for i = 0 to N − 1
    α0(i) = c0α0(i)
next i

// 각 관측 αt(i)에서의 전방계산  
for t = 1 to T − 1 # 모든 관측열에 대해서 
    ct = 0 # 해당 시간t에서의 전방변수(확률들의 합) 
    for i = 0 to N − 1 # 모든 상태에 대해서 
        αt(i) = 0  # t번째 관측에 대해서 전방계산 
        for j = 0 to N − 1 # 모든 상태에 대해서 
            αt(i) = αt(i) + αt−1(j)aji # 이전전방확률 * j -> i 상태로 전이확률의 합
        next j
        αt(i) = αt(i)bi(Ot) #t(t에서i상태)까지의 전방확률 * t시간의 i상태에서 관측된O의 확률 
        ct = ct + αt(i) # 현재 시간(t) 에서의 전방변수 (모든 상태에서의 확률의 합) 

    next i

    // αt(i) 에 대한 스케일 조정
    ct = 1/ct
    for i = 0 to N − 1
        αt(i) = ctαt(i) 
    next i
next t


3. The β-pass

// Let βT −1(i) = 1, scaled by cT −1
for i = 0 to N − 1 # 모든 상태에 대해서 
    βT −1(i) = cT −1  
next i

// β-pass
for t = T − 2 to 0 by − 1 # 끝에서 2번째 시간(관측) ~ 첫시간(관측) 까지 하나씩 감소
    for i = 0 to N − 1 # 모든 상태에 대해서 
        βt(i) = 0 # t 시간에서 i 상태인 경우에 대한 b-pass변수 
        for j = 0 to N − 1 # 모든 상태에 대해서 (상태전이가 N*N 만큼일어나니깐~i*j하는것) 
            # 현재t(i)에서의 b-pass 는 상태전이*다음시간관측확률*다음시간이후의b-pass
            βt(i) = βt(i) + aij bj (Ot+1)βt+1(j) 
        next j

        // scale βt(i) with same scale factor as αt(i)
        βt(i) = ctβt(i)
    next i
next t

# α-pass 와  β-pass 를 구했고 아래에서 그 둘을 이용해서 hidden 변수를 구한다. 
4. Compute γt(i, j) and γt(i)

for t = 0 to T − 2 # 모든 시간(관측열)에 대해서 각각 계산한다. (마지막 제외) 
    denom = 0

    for i = 0 to N − 1 # 모든 상태에 대해서
        for j = 0 to N − 1 # 모든 상태에 대해서 (N*N)
            # t에서의 전방변수(i이전)*전이확률(i*j)*다음상태(j)에서 O관측확률*후방변수(j이후) 
            denom = denom + αt(i)aij bj (Ot+1)βt+1(j)
        next j
    next i

    # demon 은 특정시간 t에서의 모든 전,후방확률들의 합 
    for i = 0 to N − 1  # 모든 상태에 대해서
        γt(i) = 0
        for j = 0 to N − 1  # 모든 상태에 대해서
            # 모든 상황(demon) 으로 나누어줌, 즉 i->j 로만 가는 확률 
            γt(i, j) = (αt(i)aij bj (Ot+1)βt+1(j))/denom 
            # γt(i) 는 결국 i 상태에서 변화되는 모든 γt(i, j) 확률들의 합이다.
            γt(i) = γt(i) + γt(i, j)  
        next j
    next i
next t

// Special case for γT −1(i)
denom = 0
for i = 0 to N − 1
    denom = denom + αT −1(i)
next i
for i = 0 to N − 1
    γT −1(i) = αT −1(i)/denom
next i


# γt(i) 는 t시간에 i상태에서 변화되는 모든 확률의 합이다. i->0, i->1, i->j ..i->N

# hidden 변수 γt(i, j) and γt(i) 를 이용하여 (A,B,π) 를 재추정한다. 
5. Re-estimate A, B and π

# π 를 재추정한다.
for i = 0 to N − 1  # 모든 상태의 개수만큼 재추정, 즉 어떤 상태로 초기화가 되겠냐?
    πi = γ0(i) # 0시간(첫번째관측)의 i상태에 있을 확률 
next i

# A (상태전이확률행렬) 재추정 (i 에서 출발하는 모든 확률로 i->특정j 확률을 나누어서 구함)
for i = 0 to N − 1  # 모든 상태의 개수만큼 재추정
    for j = 0 to N − 1  # 모든 상태의 개수만큼 재추정 (상태변이행렬은 N*N이다.)
        numer = 0
        denom = 0
        for t = 0 to T − 2  # 모든 시간(관측열)에 대해서 계산한다.
            numer = numer + γt(i, j) # 특정 시간t에서 i->j 경우의 γ 를 분자로  
            denom = denom + γt(i)  # 특정 시간t에서 i 상태인 경우의 γ 를 분모로
        next t
        aij = numer/denom  # i->j 경우 / i 에서 이동하는 모든 경우 
    next j
next i

// B(관측확률행렬) 재추정 (특정 상태에서 관측될 모든 확률 구함, N*M 이 되겠지)
for i = 0 to N − 1  # 모든 상태의 개수만큼 재추정
    for j = 0 to M − 1  #모든 관측의 개수만큼 재추정
        numer = 0
        denom = 0
        for t = 0 to T − 1  # 모든 시간(관측열)에 대해서 계산한다.
            if(Ot == j) then 
                numer = numer + γt(i) #i상태에서 관측열에서의 O가 관측될 확률 
            end if

            denom = denom + γt(i) #i상태에서 관측되는 모든것들 합 
        next t
        bi(j) = numer/denom  #i상태에서 특정O가 관측될 확률/i상태에서 아무거나 관측될 확률    
    next j
next i


# 새로 평가된 상태전이확률,관측확률,초기확률을 가지고 관측열을 로그우도로 평가한다.
6. Compute log[P(O | λ)]

logProb = 0
for i = 0 to T − 1
    logProb = logProb + log(ci)
next i
logProb = −logProb


# 평가결과를 보고 멈출지 계속 순회할지 결정한다. 
7. To iterate or not to iterate, that is the question. . .

iters = iters + 1
if (iters < maxIters and logProb > oldLogProb) then
    oldLogProb = logProb
    goto 2
else
    output λ = (π, A, B)
end if

3. 수학식으로 정리 (오일석저 '패턴인식' 에서 발췌) 

a-pass 와 b-pass  :  이전 a-pass * 현재상태로전이 * 현재상태에서 관측확률 로 이루어진다.

감마(i)   :  a-pass  와 b-pass 를 이용해 간단히 만들어진다. k-means 의 경우 E 구간에서 포인트들이 이동하는 모습들이 눈에 보이기 때문에 유추하기 쉽지만, HMM 에서 이걸 어떻게 생각했는지 참 대단하다는 생각이 든다.

감마(i,j) :  역시 a-pass  와 b-pass 를 이용해 간단히 만들어진다. 상태변이확률이 더해졌다.

new A (상태전이확률) :  감마(i) 와 감마(i,j) 로 만들어진다. 여기서는 k(i,j) 가 감마(i,j) 이다.

new B (관측확률) :  감마(i) 로만들어진다.

4. Python 코드 

# HMM의 3가지 문제 ('평가','디코딩','학습') 중 '평가'  
def evaluate(self, sequence):
    length = len(sequence)
    if length == 0:
        return 0

    prob = 0
    alpha = self._forward(sequence)
    for state in alpha[length - 1]:
        prob += alpha[length - 1][state] # 전방계산에서 마지막 시간에서 모든 상태확률의 합

    return prob

# 전방변수 구하기 α-pass
def _forward(self, sequence):  
    sequence_length = len(sequence)  #관측열의 갯수 
    if sequence_length == 0:
        return []

    alpha = [{}]    # αt(i) ,a-pass 에서 마지막 시간의 상태별 확률들의 합은 평가의 값이다. 
    
    # 초기확률 구하기 
    for state in self._states:
     # 초기상태확률 * 각 상태에서의 첫번째 관측에 대한 관측확률 
        alpha[0][state] = self.start_prob(state) * self.emit_prob(state, sequence[0])

  # 이 후의 모든 확률 구하기
    for index in xrange(1, sequence_length): # 모든 관측열(시간T) 에 대해서 분석 
        alpha.append({})
        for state_to in self._states: # 모든 상태로 도착 
            prob = 0
            for state_from in self._states: # 모든 상태에서 시작 
          # 이전의 a-pass 과 to상태로의 상태전이확률의 합이 prob 이 된다.
                prob += alpha[index - 1][state_from] * self.trans_prob(state_from, state_to)
        
       # 현재 상태에 대한 a-pass 값은 prob * to상태에서의 관측확률
            alpha[index][state_to] = prob * self.emit_prob(state_to, sequence[index])

    return alpha

# 후방변수 구하기 β-pass 
def _backward(self, sequence):
    sequence_length = len(sequence)
    if sequence_length == 0:
        return []

    beta = [{}]
    for state in self._states:
        beta[0][state] = 1

    for index in xrange(sequence_length - 1, 0, -1):
        beta.insert(0, {})
        for state_from in self._states:
            prob = 0
            for state_to in self._states:
                prob += beta[1][state_to] * \
                        self.trans_prob(state_from, state_to) * \
                        self.emit_prob(state_to, sequence[index])
            beta[0][state_from] = prob

    return beta

"""
    주어진 `sequence` 로 적절한 상태전이행렬과 관측행렬을 구한다.
    `smoothing`은 기본적으로 0이며 관련내용은은 아래 링크에서 확인하라.
    `additive smoothing <http://en.wikipedia.org/wiki/Additive_smoothing>`
"""
# 학습 시작 
def learn(self, sequence, smoothing=0):

    length = len(sequence)
    alpha = self._forward(sequence)  # 전방 a-pass 구하기
    beta = self._backward(sequence) # 후방 b-pass 구하기 

  # γt(i) 를 구한다. 간단히 말해 (a-pass * b-pass ) / sum 
    gamma = []
    for index in xrange(length):   # 모든 관측열(시간T) 에 대해서 계산해준다. 
        prob_sum = 0
        gamma.append({})
        for state in self._states:  # 모든 관측열(시간T)의  모든 상태에 대해서 계산해준다.
       # 특정 시간/상태에서의 a-pass * 특정 시간/상태에서의 b-pass
            prob = alpha[index][state] * beta[index][state]
            gamma[index][state] = prob  
            prob_sum += prob

        if prob_sum == 0:
            continue

        for state in self._states:  # 모든 상태에 대해서 
            gamma[index][state] /= prob_sum  # 특정시간/상태에 대한 전체로 나눈 감마 확률을 구한다.

    xi = []
    for index in xrange(length - 1):
        prob_sum = 0
        xi.append({})
        for state_from in self._states:
            xi[index][state_from] = {}
            for state_to in self._states:
                prob = alpha[index][state_from] * beta[index + 1][state_to] * \
                       self.trans_prob(state_from, state_to) * \
                       self.emit_prob(state_to, sequence[index + 1])
                xi[index][state_from][state_to] = prob
                prob_sum += prob

        if prob_sum == 0:
            continue

        for state_from in self._states:
            for state_to in self._states:
                xi[index][state_from][state_to] /= prob_sum

    states_number = len(self._states)
    symbols_number = len(self._symbols)
    for state in self._states:
        # update start probability
        self._start_prob[state] = \
            (smoothing + gamma[0][state]) / (1 + states_number * smoothing)

        # update transition probability
        gamma_sum = 0
        for index in xrange(length - 1):
            gamma_sum += gamma[index][state]

        if gamma_sum > 0:
            denominator = gamma_sum + states_number * smoothing
            for state_to in self._states:
                xi_sum = 0
                for index in xrange(length - 1):
                    xi_sum += xi[index][state][state_to]
                self._trans_prob[state][state_to] = (smoothing + xi_sum) / denominator
        else:
            for state_to in self._states:
                self._trans_prob[state][state_to] = 0

        # update emission probability
        gamma_sum += gamma[length - 1][state]
        emit_gamma_sum = {}
        for symbol in self._symbols:
            emit_gamma_sum[symbol] = 0

        for index in xrange(length):
            emit_gamma_sum[sequence[index]] += gamma[index][state]

        if gamma_sum > 0:
            denominator = gamma_sum + symbols_number * smoothing
            for symbol in self._symbols:
                self._emit_prob[state][symbol] = \
                    (smoothing + emit_gamma_sum[symbol]) / denominator
        else:
            for symbol in self._symbols:
                self._emit_prob[state][symbol] = 0

전체소스: https://github.com/jason2506/PythonHMM/blob/master/hmm.py 

레퍼런스:

https://www.cs.sjsu.edu/~stamp/RUA/HMM.pdf
http://article.sciencepublishinggroup.com/html/10.11648.j.acm.s.2017060401.12.html#paper-content-1
https://github.com/nilmtk/nilmtkhttps://github.com/jason2506/PythonHMM/blob/master/hmm.py 

 * 참고로 "밑바닥부터 배우는 데이터과학" 서적은 numpy,scikit-learn 등의 외부라이브러리를 활용은 배제되었습니다.

def least_squares_fit(x,y):
    beta = correlation(x, y) * standard_deviation(y) / standard_deviation(x)
    alpha = mean(y) - beta * mean(x)
    return alpha, beta

def predict(alpha, beta, x_i):
    return beta * x_i + alpha

def error(alpha, beta, x_i, y_i):
    return y_i - predict(alpha, beta, x_i)

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

def squared_error_gradient(x_i, y_i, theta):
    alpha, beta = theta
    return [-2 * error(alpha, beta, x_i, y_i),       # alpha partial derivative
            -2 * error(alpha, beta, x_i, y_i) * x_i] # beta partial derivative

random.seed(0)
theta = [random.random(), random.random()]
alpha, beta = minimize_stochastic(squared_error, 
                                  squared_error_gradient,
                                  num_friends_good,
                                  daily_minutes_good, 
                                  theta,
                                  0.0001)

def minimize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):

    data = zip(x, y)
    theta = theta_0                             # initial guess
    alpha = alpha_0                             # initial step size
    min_theta, min_value = None, float("inf")   # the minimum so far
    iterations_with_no_improvement = 0

    # if we ever go 100 iterations with no improvement, stop
    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:
            # if we've found a new minimum, remember it
            # and go back to the original step size
            min_theta, min_value = theta, value
            iterations_with_no_improvement = 0
            alpha = alpha_0
        else:
            # otherwise we're not improving, so try shrinking the step size
            iterations_with_no_improvement += 1
            alpha *= 0.9

        # and take a gradient step for each of the data points
        for x_i, y_i in in_random_order(data):
            gradient_i = gradient_fn(x_i, y_i, theta)
            theta = vector_subtract(theta, scalar_multiply(alpha, gradient_i))

    return min_theta

레퍼런스 :

밑바닥부터 시작하는 데이터 과학

만화로 쉽게 배우는 회귀분석  

요즘 데이터분석과 관련하여  텐서플로우와 스파크(ML)등의 머신러닝 솔루션들이 굉장히 유행하고 있습니다. 물론 저것들이 삶을 편안하게 만들어주기도 하지만 대부분의 데이터 분석은 저런 거창한 것 말고 평균,편차 같은 기본적인 개념으로 부터 시작되고 있으며 이러한 개념을 조금씩 변경해가며 더 의미있는 가치를 찾기 위해 빠르게 적용해보는 과정을 거쳐야하는데  그러기 위해서는 

1. 직접 코딩해서 기본적인 데이터분석 유틸리티 함수들을 만들어봐야한다. 

2. SQL문을 잘 다루어야한다. 

3. 엑셀을 잘 다루어야 한다. 

이 3가지는 기본이라고 저는 생각합니다. 소규모 데이터를 가지고 이리저리 반복해서 돌려보는 과정은 매우 중요하며 이런 기본적인 것들도 못하면서 하둡,텐서플로우나 깔짝대고, 데이터 분석 한다 라고 칭할 수는 없겠지요. (이것은 논쟁의 여지가 있습니다.) 그래서 이것들 중 1번에 대하여  "밑바닥부터 시작하는 데이터과학" 등의 좋은책의 내용을 통해서 살펴보는 시간을 갖도록 하겠습니다. 

통계,확률,패턴인식 분야의 내용들의 수식은 외우기가 힘듭니다. 잘은 모르지만 외워도 금방 까먹어지는게 통계관련 공식인거 같습니다. (예를들어 정규분포,확률밀도라는 개념은 매우 쉽게 수긍이 가지만 , 그 식을 외우는건 좀 .;;) 또한 체득하고 나면 당연한거 아냐? 라고 느껴지는 알고리즘(수식) 인데 글로 설명하면 매우 산만해지는 학문인거 같습니다. 하지만  우리들은 소프트웨어 개발자이기 때문에 수많은 기술변화도 따라가야하는 운명에 있는데 쓸때 없이 공식 및 수학을 외우고 있을 수는 없지 않겠습니까? 

제가 초등학교 5학년때 만(10,000) 단위 암산을 했었는데요. 주산,암산을 배워보신분은 아시겠지만 그것이 가능한것은 기억의 매개로 주판을 이용한 것 입니다. 마찬가지로 우리 개발자들은 그 기억의 매개체로 코드를 이용하면 개념과 함정(예를들어 K-Means 를 통해 군집화하면 길이차를 가지고 구분짓는 것이기 때문에 문제가 생기는 도메인 또한 많다) 에 대한 이해가 더 오래갈 것이라 보며 덤으로 가져다 사용도 할 수 있을것입니다.

* 그 상상의 매개체로써의 언어로 "파이썬" 을 선택했으며 정말 좋은 언어라고 생각합니다.

 * 참고로 "밑바닥부터 배우는 데이터과학" 서적은 numpy,scikit-learn 등의 외부라이브러리를 활용은 배제되었습니다.

단점은 위의 그림처럼 엄한데 가서 최소값을 찾을 수 도 있다는 것이다. 진정한 최소 값은 다른 곳에 있을 수도 있다.

코드로 말해요.

함수하나를 만들어보자. 실수 백터를 입력하면  요소의 제곱의 합을  리턴해 주는 함수이다. (비용함수라고 치자)

def sum_of_squares(v):
    return sum(v_i ** 2 for v_i in v)

위에서 v 는 아래 설명식에서의 t (1~N) 에서의 X 값과 같다) 물론 동일식은 아니다.

미분값을 구해보자.  

f 라는 함수에 대해 x 위치에서의 미분값 즉 기울기를 리턴한다. (위의 식에서의 x가 아니다. 헥깔리지 마시라) h 는 수학에서는 극한으로 가야겠지만 프로그래밍에서는 근사로 충분하기 때문에 적절히 작은수(0.00001 정도?) 를 넣어준다. 

이 미분값의 방향에 따라서 조금씩 스텝을 바꿔서 최소값으로 나아가게 된다.
(여기서 x 는 위의 해설에서 m 과 b 와 같다, 즉 여기서는 x 를 조절해서 최소의 기울기를 찾아나간다) 

def difference_quotient(f, x, h):
    return (f(x + h) - f(x)) / h

여러 변수중 하나의 변화율만 관찰하는 편도함수를 구해보자.  

def partial_difference_quotient(f, v, i, h):
    """ 함수 f의 i번째 편도함수가 v에서 가지는 값 """

    w = [v_j + (h if j == i else 0)      # h를 v의 i번째 변수에만 더해주자.
         for j, v_j in enumerate(v)]      # 즉 i 번째 변수만 변화할 경우 

    return (f(w) - f(v)) / h

i 번째 변수에 대한 편 변화율을 계산한다. 즉 i번째 편도함수는 , i 번째 변수를 제외한 다른 모든 입력변수를 고정시켜서 계산한다. 

def estimate_gradient(f, v, h=0.00001):
    return [partial_difference_quotient(f, v, i, h)
            for i, _ in enumerate(v)]

모든 변수에 대한 편 변화율을 구한다. (여기서 v는 해설 예에서 m,b 이다.) 

Gradient 적용하기 

아래의 sum_of_squares_gradient 는 이미 정의된 도함수이다.  (x^2  ->  2x 로 미분)

# 8.3. 경사 적용하기

def step(v, direction, step_size):
    """ move step_size in the direction from v"""

    return [v_i + step_size * direction_i
            for v_i, direction_i in zip(v, direction)]


def sum_of_squares_gradient(v):
    return [2 * v_i for v_i in v]

#임의의 시작점을 선택 
v = [random.randint(-10,10) for i in range(3)]
tolerance = 0.0000001

while True:
    #print v, sum_of_squares(v)
    gradient = sum_of_squares_gradient(v) # v 에서의 경사기울기(gradient)를 구한다.
    next_v = step(v, gradient, -0.01)     # 그 기울기의 반대방향으로 이동
    if distance(next_v, v) < tolerance:   # 만약 기준내로 이동하였다면 멈춤. 
        break
    v = next_v                            # 그렇지 않다면 다시 이동.

적절한 이동 거리 정하기

너무 많이 이동하면 발산될 위험이 있고, 너무 적게 이동하면 계산시간이 너무 오래 걸린다.

적절한 거리를 미리 정해두고 가장 합리적인 값을 선택한다.

step_sizes = [100, 10, 1, 0.1, 0.01, 0.001, 0.0001, 0.00001]

def safe(f):
    """define a new function that wraps f and return it"""
    def safe_f(*args, **kwargs):
        try:
            return f(*args, **kwargs)
        except:
            return float('inf')         # this means "infinity" in Python
    return safe_f

종합하기

def minimize_batch(target_fn, gradient_fn, theta_0, tolerance=0.000001):
    """use gradient descent to find theta that minimizes target function"""

    step_sizes = [100, 10, 1, 0.1, 0.01, 0.001, 0.0001, 0.00001]

    theta = theta_0                           # set theta to initial value
    target_fn = safe(target_fn)               # safe version of target_fn
    value = target_fn(theta)                  # value we're minimizing

    while True:
        gradient = gradient_fn(theta)
        next_thetas = [step(theta, gradient, -step_size)
                       for step_size in step_sizes]

        # choose the one that minimizes the error function
        next_theta = min(next_thetas, key=target_fn)
        next_value = target_fn(next_theta)

        # stop if we're "converging"
        if abs(value - next_value) < tolerance:
            return theta
        else:
            theta, value = next_theta, next_value

SGD  (Stochastic gradient descent) 

앞의 minimize_batch 를 이용하면 반복문을 돌 때마다 데이터 전체에 대해 gradient 값을 계산해야 하기 때문에 계산시간이 오래걸린다. 그런데 대부분의 오류 함수는 더 할 수 있는 속성을 갖고 있다. 즉 데이터 전체에 대한 오류값이 각각 데이터 포인트에 대한 오류값의 합과 같다. 이럴 때는 한 번 반복문을 돌 때마다 데이터 포인트 한 개에 대한 gradient 를 계산해야하는 SGD 를 사용할 수 있다.

def in_random_order(data):
    """generator that returns the elements of data in random order"""
    indexes = [i for i, _ in enumerate(data)]  # create a list of indexes
    random.shuffle(indexes)                    # shuffle them
    for i in indexes:                          # return the data in that order
        yield data[i]


def minimize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):

    data = zip(x, y)
    theta = theta_0                             # initial guess
    alpha = alpha_0                             # initial step size
    min_theta, min_value = None, float("inf")   # the minimum so far
    iterations_with_no_improvement = 0

    # if we ever go 100 iterations with no improvement, stop
    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:
            # if we've found a new minimum, remember it
            # and go back to the original step size
            min_theta, min_value = theta, value
            iterations_with_no_improvement = 0
            alpha = alpha_0
        else:
            # otherwise we're not improving, so try shrinking the step size
            iterations_with_no_improvement += 1
            alpha *= 0.9

        # and take a gradient step for each of the data points
        for x_i, y_i in in_random_order(data):
            gradient_i = gradient_fn(x_i, y_i, theta)
            theta = vector_subtract(theta, scalar_multiply(alpha, gradient_i))

    return min_theta

레퍼런스 :

밑바닥부터 시작하는 데이터 과학

 * 참고로 "밑바닥부터 배우는 데이터과학" 서적은 numpy,scikit-learn 등의 외부라이브러리를 활용은 배제되었습니다.

데이터 다루기

대부분의 언어들이 함수형을 차용하고 있으며 , 파이썬도 빠질리는 없다. 파이썬의 함수형 도구들을 사용하여 데이 터 분석을 해보자.

주요 도구: List comprehension 

함수형 (고차함수) 의 특징인 lambda , filter, map, reduce 알고리즘을 리스트로 간단히 구현하는 방법입니다.

예를 들면  

a = [1,2,3,4,5] 중에서 홀수만 골라낸다고 할때 

람다식을 이용하면 아래와 같이 됩니다.

filter (lambda x : x % 2 , a)   

근데 이것을 List comprehension 으로 표현하면

 [x for x in a if x % 2] 

이렇게 됩니다. 즉 a 라는 컬렉션에서 하나씩 값(x) 을 가져와서 , 만약 짝수면 리스트의 요소로 추가합니다.

파이썬에서는 좀 황당한게 List comprehension 이 더 사용성 (가독성) 이 좋아서 오래전 부터 고차함수들을 제꼈(?)습니다. @@

코드로 말해요. 

@ 아래와 같은 주식 데이터가 있다고 하자.  / 종가 / 날짜 / 심볼로 이루어져있다.

data = [        {'closing_price' : 102.06,
'date' : datetime.datetime(2014, 8, 29, 0, 0),
'symbol' : 'AAPL'}, 
... 
...

* 파일에 저장된 모습 (아래 다운로드 받을 수 있다)

symbol  date   closing_price
AAPL   2015-01-23 112.98
AAPL   2015-01-22 112.4
AAPL   2015-01-21 109.55
AAPL   2015-01-20 108.72
AAPL   2015-01-16 105.99
AAPL   2015-01-15 106.82
AAPL   2015-01-14 109.8
..

stocks.txt

@ 심볼이 'AAPL' 인 것들 중에서  가장 높은 가격 산출

max_aapl_price = max(row["closing_price"]
                              for row in data if row["symbol"] == "AAPL")

print "max aapl price", max_aapl_price

data 컬렉션에서 값(row) 을 하나 가져와서  심볼이 "AAPL" 인것들 중에서 값이 가장 높은것

@ 심볼별로 데이터 그룹화

# group rows by symbol
by_symbol = defaultdict(list)

for row in data:
    by_symbol[row["symbol"]].append(row)

심볼이 같은것 끼리 묶어 줍니다. 즉 { "AAPL" : [ ... , ... ]  , "BBPL" : [ ..., ... ]  .......} 이런식으로요~

defaultdict 은 파이썬에서 제공하는 컬렉션 중에 하나로써 , 만약 키 값이 없더라도 새로 생성해서 append 하게 해 줍니다. 일반 딕셔너리로 저렇게 코딩하면 에러죠. 

@ 그룹별 가장 높은 가격 

max_price_by_symbol = {symbol: max(row["closing_price"] for row in grouped_rows)
                                 for symbol, grouped_rows in by_symbol.iteritems()}

위에서 심볼별 리스트들을 가져와서 각각에 대한 가장 높은 값을 구해줍니다.

결과는 {'AAPL': 119.0, 'FB': 38.23, 'MSFT': 22.78}

@ 특정필드를 리스트로 추출  

def picker(field_name):
    return lambda row: row[field_name]

def pluck(field_name, rows):
    return map(picker(field_name), rows)

picker 는 특정 필드의 값을 리턴해주는 함수를 반환합니다.

pluck 는 rows 들 중 특정 필드의 값 만으로 리스트를 만들어 줍니다.  

즉 field_name 에 "closing_price" 를 입력해주면 , 그 값들로만 이루어진 리스트를 리턴함. 

@ 데이터를 그룹화시켜서 모으고 , value_transform 공식에 의해서 변경시킨다. 

def group_by(grouper, rows, value_transform=None):
    grouped = defaultdict(list)
    for row in rows:
        grouped[grouper(row)].append(row)
    if value_transform is None:
        return grouped
    else:
        return {key: value_transform(rows)
                for key, rows in grouped.iteritems()}

심볼이라든지 어떤 특정 것을 구분하여 그룹화시키고 , 각 그룹을 value_transform 에 의해 계산합니다.
value_transform  가 그룹의 max 값을 찾는 것이라면 {"AAPL" : 1030 , "BBPL" : 2000 ..} 뭐 이런식으로 결과가 나오겠지요.

@ 그룹된 로우를 date 에 의해 정렬한 후에 이전날과 비교한 오늘의 종가의 증감비율을 추가한다.  

위에 언급한 value_tansform 로써 사용된다. 

def percent_price_change(yesterday, today):
    return today["closing_price"] / yesterday["closing_price"] - 1

def day_over_day_changes(grouped_rows):
    # sort the rows by date
    ordered = sorted(grouped_rows, key=picker("date"))
    # zip with an offset to get pairs of consecutive days
    return [{"symbol": today["symbol"],
             "date": today["date"],
             "change": percent_price_change(yesterday, today)}
            for yesterday, today in zip(ordered, ordered[1:])]

@ 데이터를 "심볼" 별로 그룹핑하고, 이전날과 비교한 오늘의 종가의 증감비율을 계산하 day_over_say_changes 를  value_tansform 로 넘겨준다. change 키에 그 값들을 넣어준다.

changes_by_symbol = group_by(picker("symbol"), data, day_over_day_changes)

@ 모든 그룹들의  "Change"  의 값들중 max,min 구해서 하나의 거대한 리스트에  담는다.

all_changes = [change
               for changes in changes_by_symbol.values()
               for change in changes]

print "max change", max(all_changes, key=picker("change"))
print "min change", min(all_changes, key=picker("change"))

data = [[1, 20, 2],
        [1, 30, 3],
        [1, 40, 4]]
def shape(A):
    num_rows = len(A)
    num_cols = len(A[0]) if A else 0
    return num_rows, num_cols

def scale(data_matrix):
    num_rows, num_cols = shape(data_matrix)
    means = [mean(get_column(data_matrix,j))
             for j in range(num_cols)]
    stdevs = [standard_deviation(get_column(data_matrix,j))
              for j in range(num_cols)]
    return means, stdevs

def rescale(data_matrix):
    means, stdevs = scale(data_matrix)

    def rescaled(i, j): 
        if stdevs[j] > 0:
            return (data_matrix[i][j] - means[j]) / stdevs[j]
        else:
            return data_matrix[i][j]

    num_rows, num_cols = shape(data_matrix)
    return make_matrix(num_rows, num_cols, rescaled)

X = [
    [20.9666776351559, -13.1138080189357],
    [22.7719907680008, -19.8890894944696],
    [25.6687103160153, -11.9956004517219],
    [18.0019794950564, -18.1989191165133],
     ....
  
    [15.6563953929008, -17.2196961218915],
    [25.2049825789225, -14.1592086208169]
]

   def shape(A):
      num_rows = len(A)
      num_cols = len(A[0]) if A else 0
      return num_rows, num_cols
 def magnitude(v):
    return math.sqrt(sum_of_squares(v))

   def de_mean_matrix(A):
        """A 행렬의 모든 값에 해당 컬럼의 평균을 뺀 행렬을 리턴합니다"
        nr, nc = shape(A)
        column_means, _ = scale(A)
        return make_matrix(nr, nc, lambda i, j: A[i][j] - column_means[j])

def negate(f):
    """-f(x) 로 리턴한다. """
    return lambda *args, **kwargs: -f(*args, **kwargs)


def negate_all(f):
    """모든 리턴을 -y 로 한다."""
    return lambda *args, **kwargs: [-y for y in f(*args, **kwargs)
# 8.4. 적절한 스텝 크기 정하기

def safe(f):
    """f 를 감싸고 it 를 리턴하는 새 함수를 정의한다."""
    def safe_f(*args, **kwargs):
        try:
            return f(*args, **kwargs)
        except:
            return float('inf')         # this means "infinity" in Python
    return safe_f

def minimize_batch(target_fn, gradient_fn, theta_0, tolerance=0.000001):
    """ 경사하강법을 이용하여 타겟 함수를 최소화할 세타를 찾는다."""

    step_sizes = [100, 10, 1, 0.1, 0.01, 0.001, 0.0001, 0.00001]

    theta = theta_0                           # 세타 초기값 
    target_fn = safe(target_fn)               # target_fn 의 safe 버전 
    value = target_fn(theta)                  # 최소화할 값 

    while True:
        gradient = gradient_fn(theta)
        next_thetas = [step(theta, gradient, -step_size)
                       for step_size in step_sizes]

        # choose the one that minimizes the error function
        next_theta = min(next_thetas, key=target_fn)
        next_value = target_fn(next_theta)

        # stop if we're "converging"
        if abs(value - next_value) < tolerance:
            return theta
        else:
            theta, value = next_theta, next_value

def maximize_batch(target_fn, gradient_fn, theta_0, tolerance=0.000001):
    return minimize_batch(negate(target_fn),
                          negate_all(gradient_fn),
                          theta_0,
                          tolerance)
# partial 은 functools 에서 제공하는 함수로 
# https://www.pydanny.com/python-partials-are-fun.html 참고 
def first_principal_component(X):
    guess = [1 for _ in X[0]]
    unscaled_maximizer = maximize_batch(
        partial(directional_variance, X),  # is now a function of w
        partial(directional_variance_gradient, X),  # is now a function of w
        guess)
    return direction(unscaled_maximizer)

def principal_component_analysis(X, num_components):
    components = []
    for _ in range(num_components):
        component = first_principal_component(X)
        components.append(component)
        X = remove_projection(X, component)
        
    return components

Y = de_mean_matrix(X)
components = principal_component_analysis(Y, 2)
print "principal components", components
print "first point", Y[0]
print "first point transformed", transform_vector(Y[0], components)

레퍼런스:

밑바닥부터 배우는 데이터과학

요즘 데이터분석과 관련하여  텐서플로우와 스파크(ML)등의 머신러닝 솔루션들이 굉장히 유행하고 있습니다. 물론 저것들이 삶을 편안하게 만들어주기도 하지만 대부분의 데이터 분석은 저런 거창한 것 말고 평균,편차 같은 기본적인 개념으로 부터 시작되고 있으며 이러한 개념을 조금씩 변경해가며 더 의미있는 가치를 찾기 위해 빠르게 적용해보는 과정을 거쳐야하는데  그러기 위해서는 

1. 직접 코딩해서 기본적인 데이터분석 유틸리티 함수들을 만들어봐야한다. 

2. SQL문을 잘 다루어야한다. 

3. 엑셀을 잘 다루어야 한다. 

이 3가지는 기본이라고 저는 생각합니다. 소규모 데이터를 가지고 이리저리 반복해서 돌려보는 과정은 매우 중요하며 이런 기본적인 것들도 못하면서 하둡,텐서플로우나 깔짝대고, 데이터 분석 한다 라고 칭할 수는 없겠지요. (이것은 논쟁의 여지가 있습니다.) 그래서 이것들 중 1번에 대하여  "밑바닥부터 시작하는 데이터과학" 등의 좋은책의 내용을 통해서 살펴보는 시간을 갖도록 하겠습니다. 

통계,확률,패턴인식 분야의 내용들의 수식은 외우기가 힘듭니다. 외워도 금방 까먹어지는게 통계관련 공식인거 같습니다. (예를들어 정규분포,확률밀도라는 개념은 매우 쉽게 수긍이 가지만 , 그 식을 외우는건 좀 .;;) 또한 체득하고 나면 당연한거 아냐? 라고 느껴지는 알고리즘(수식) 인데 글로 설명하면 매우 산만해지는거 학문인거 같습니다. 하지만  우리들은 소프트웨어 개발자이기 때문에 수많은 기술변화도 따라가야하는 운명에 있는데 쓸때 없이 공식 및 수학을 외우고 있을 수는 없지 않겠습니까? 

제가 초등학교 5학년때 만(10,000) 단위 암산을 했었는데요. 주산,암산을 배워보신분은 아시겠지만 그것이 가능한것은 기억의 매개로 주판을 이용한 것 입니다. 마찬가지로 우리 개발자들은 그 기억의 매개체로 코드를 이용하면 개념과 함정(예를들어 K-Means 를 통해 군집화하면 길이차를 가지고 구분짓는 것이기 때문에 문제가 생기는 도메인 또한 많다) 에 대한 이해가 더 오래갈 것이라 보며 덤으로 가져다 사용도 할 수 있을것입니다.

* 그 상상의 매개체로써의 언어로 "파이썬" 을 선택했으며 정말 좋은 언어라고 생각합니다.

 * 참고로 "밑바닥부터 배우는 데이터과학" 서적은 numpy,scikit-learn 등의 외부라이브러리를 활용은 배제되었습니다.

연관 되는 이벤트들이란?

이마트에 가서 장을 보는 것으로 이야기를 풀어나가 본다.
어느 날이 좋은 목요일엔  ("수박"-"기저귀"-"맥주") 를 사고 어느날 안좋은 목요일에도  ("기저귀"-"맥주") 를 샀다고하자. 그런 날이 자주 있었다고 하면 목요일과 "기저귀" 와 "맥주" 는 어떤 연관관계가 있다고 할 수 있다. 나 뿐만 아니라 수많은 남자들이 목요일날 저렇게 함께 사는 확률이 높다는것을 알아차린다면 업주는 무엇을 해야할까? 그렇다 "기저귀" 와 "맥주" 코너를 비슷한 곳에 위치시킴으로써 매출을 늘릴 수 있을 것이다.

이렇게 어떤 데이터 집합속에서 각 데이터간의 관계를 살펴서 연관 되는 분석/규칙을 찾는 알고리즘중에 하나가 Apriori 알고리즘이다. 먼저 관련 용어를 아래 표를 가지고 설명하겠다.  매우 간단하다.

지지도 

위의 표로 보면 [두유] 의 지지도는 4/5 인데 이유는 전체 집합군에서 두유가 포함된 집합의 수이다.

전체에서 아이템 집합 (두유) 이 포함된 데이터 집합의 비율로 정의한다. 

[두유,기저귀] 의 지지도는 ? 그렇다. 3/5 이다.  

두유와 기저귀를 모두 포함한 데이터 집합은 3개가 있다.

신뢰도 

[기저귀]를 샀을때 [맥주]도 같이 사는 확률에 관한 이야기이다.

즉 [기저귀] ->[맥주] 로 그 둘간의 연관 규칙을 나타내는데,  

[기저귀] 를 산 모든날이 100일이라고 하면 , [기저귀] 와 [맥주] 를 함께 산 날이 90일이라고 칠때

그 신뢰도는 9/10이다. 즉 신뢰도란 지지도({기저귀,맥주})  / 지지도({기저귀}) 가 된다.

Apriori 알고리즘이란?

위의 그림은 A 가 발생하고 나서 (A,B) (A,C) (A,D) 등이 일어날 수 있고 그 후에 (A,B,C) 가 일어날수 있음에 대한 단순한 트리이다. 이때 AB가 일어날 확률이 적었다면 , AB에서 파생되는 것들도 적을것임은 너무 자명하다. 따라서그것들에 대한 계산을 삭제해가면서 빈발(Frequent) 관계를 찾는 알고리즘이라고 보면된다. 알고리즘의 성능을 높게하기 위해서 말이다. 

소스로 말해요.

이제 알고리즘을 만들기 위해 필요한 함수들을 공부해보도록 하자. 파이썬 공부를 겸하면서~1석2조!!

[[1,3,4], [2,3,5], [1,2,3,5], [2,5]] 이런 데이타가 있다고 하자.

1. createC1 

데이터셋들을 중복되지 않게 유일한 수들로 나열시키는 함수이다.
[1,1,1,2,3,2,2,3,3,4,5,3,4] 의 결과는 [1,2,3,4,5] 로 나올것이다.

def createC1(dataSet):
    C1=[]
    for transaction in dataSet:
        for item in transaction:
            if not [item] in C1:
                C1.append([item])

    C1.sort()
    return map(frozenset, C1)

[1,2,3,4,5] 이렇게 나온결과를 map 을 사용하여 변경해주고 있다.
map, reduce, filter, folder 과 같은 것들은 함수형 파라다임에서 기본적인 함수들이며, 많은 언어에서 지원하는 추세이다. map 의 경우 두번째 파라미터를 , 첫번째 파라미터이 함수를 적용해서 리스트를 뽑아낸다.

결국 [frozenset([1]), frozenset([2]) .....] 이런식의 최종 결과가 도출된다. 

* frozenset 은 한번 정의하면 이후 add 같은 변경이 불가능한 집합을 의미 하며 나중에 딕셔너리의 키로 사용할 예정이기때문에 굳이 frozenset 을 사용하였다. 

2. scanD 

C1 은 [1,2,3,4,5] 같은 유일한 데이터의 리스트이고 (정확히는 [frozenset([1]), frozenset([2]) .....])
D 는 [set([1,3,4]), set([2,3,5]) .... 이다.

C1 = createC1(dataset)

D = map(set,dataset)  # distinct 수를 뽑아냄. (1,3,3) 이면 (1,3) 

# 전체 그룹중에 Ck 가 있는것.
# Ck 가 1 이라면 1을 포함한 그룹들을 찾음 (지지도 생성) 
# Ck 가 [1,3] 이라면 [1,3] 둘다 포함한 그룹들의 지지도 찾음 
def scanD(D, Ck, minSupport):
  
    # 요소 n 은 몇개의 그룹에 포함되어 있는지 계산한다. (1은 2개) 
    ssCnt = {}
    for tid in D:          # tid 에는 set([1,3,4]) 등이 담긴다.
        for can in Ck:     # can 에는 frozenset([1]) 등이 담긴다.
            if can.issubset(tid):  # 1 이 [1,3,4] 에 있으면 
                if not ssCnt.has_key(can): ssCnt[can] = 1 
                else: ssCnt[can] += 1  # 키 : 1  값 : 1을 가진 그룹의 수 

    
    # 지지도가 0.5보다 높은것들의 리스트를 구함. 
    # 1의 경우 2군데 포함되었으니 2/4 = 0.5 
    # 2의 경우 3/4 , 3의 경우 3/4 , 4의 경우 1/4, 5의 경우 3/4 
    numItems = float(len(D)) # 그룹 갯수
    retList = []
    supportData = {}
    for key , value in ssCnt.iteritems():
        support = value / numItems 
        if support >= minSupport:
            retList.insert(0,key)
        supportData[key] = support

    return retList, supportData

L1 = scanD(D, C1, 0.5) 

는 [frozenset([1]), frozenset([3]), frozenset([2]), frozenset([5])]  이런 결과가 나온다.

지지도가 0.5 이하인 4 빼고는 다 나오는셈~ 

이마트로 설명해보면 사람들이 장을 보는 아이템들이
도깨비는 [상추,기저귀,맥주,삼겹살] 를 사고
은탁이는 [상추,삼겹살] 
저승사자는 [상추]
써니는 [맥주,삼겹살] 을 샀다면  

50%이상 선택되는 아이템을 찾는 로직이며, 결과는 기저귀 빼고 나머지가 선택될 것이다.
[상추,삼겹살] 짝 또한 지지도가 50%이상 선택되어졌다. 

apriori 알고리즘 (빈발아이템 집합찾기)

# [1,3,2,5] 를 
# [1,3] [2,5] [2,3] [3,5] 를 
# [2,3,5] 이런식으로 만드는 함수 
def aprioriGen(Lk, k):
    retList = []
    lenLk = len(Lk)
    for i in range(lenLk):
        for j in range(i+1, lenLk):
            L1 = list(Lk[i])[:k-2]; L2 = list(Lk[j])[:k-2]
            L1.sort(); L2.sort()
            if L1 == L2:
                retList.append(Lk[i] | Lk[j])

    return retList

# 특정 지지도 이상의 값들의 쌍을 찾음 
def apriori(dataset, minSupport = 0.5):
    C1 = createC1(dataset)
    D = map(set, dataset)
    L1 , supportData = scanD(D,C1,minSupport)
    L = [L1]
    k=2
    while (len(L[k-2]) > 0):
        Ck = aprioriGen(L[k-2],k)
        Lk,supK = scanD(D,Ck,minSupport) # 후보그룹을 모두 찾는다.
        supportData.update(supK)
        L.append(Lk) #이게 핵심!특정 지지도 이상의 그룹들만 L에 담는다.즉 가지치기
        k += 1
    return L, supportData


if __name__ == "__main__":
   print "apriori 알고리즘"

   dataset = loadDataSet()

   L, suppData = apriori(dataset)
   print "L:"  + str(L)
   print "........................."
   print "suppData:" + str(suppData)

코드에서 retList.append(Lk[i] | Lk[j]) 는 만약 frozenset([2,3]) 과 frozenset([2,5]) 가 있다면 [2,3,5] 로 만들어주는 코드이다.  예를들어 자세히 살펴보자.

Lk = [frozenset([1,2]), frozenset([1,3])]

Lk[0] | Lk[1]  는 frozenset([1,2,3])   # 합집합

Lk[0] - Lk[1]  는 frozenset([2])       # 차집합 

Lk[0] & Lk[1]  는 frozenset([1])       # 교집합 

numbers1 = {1, 3, 5, 7}
numbers2 = {1, 3}

# # Is subset.
if numbers2.issubset(numbers1):
    print("Is a subset")

# # Is superset.
if numbers1.issuperset(numbers2):
    print("Is a superset")

# Intersection of the two sets.
   print(numbers1.intersection(numbers2))

지지도 0.5 이상의 가장 빈번한 조합들 

L:[[frozenset([1]), frozenset([3]), frozenset([2]), frozenset([5])], [frozenset([1, 3]), frozenset([2, 5]), frozenset([2, 3]), frozenset([3, 5])], [frozenset([2, 3, 5])], []]

조합들의 지지율 상세  

suppData:{frozenset([5]): 0.75, frozenset([3]): 0.75, frozenset([2, 3, 5]): 0.5, frozenset([1, 2]): 0.25, frozenset([1, 5]): 0.25, frozenset([3, 5]): 0.5, frozenset([4]): 0.25, frozenset([2, 3]): 0.5, frozenset([2, 5]): 0.75, frozenset([1]): 0.5, frozenset([1, 3]): 0.5, frozenset([2]): 0.75}

def generateRules(L, supportData, minConf=0.7):
    bigRuleList = []
    for i in range(1, len(L)):
        for freqSet in L[i]:
            H1 = [frozenset([item]) for item in freqSet]
            if i>1:
                rulesFromConseq(freqSet, H1, supportData, bigRuleList, minConf)
            else:
                calcConf(freqSet,H1,supportData, bigRuleList, minConf)

    return bigRuleList

def calcConf(freqSet, H, supportData, br1, minConf=0.7):
    prunedH = []
    for conseq in H:
        conf = supportData[freqSet] / supportData[freqSet-conseq]
        if conf >= minConf:
            print freqSet-conseq, '-->', conseq, 'conf:', conf
            br1.append((freqSet-conseq, conseq, conf))
            prunedH.append(conseq)
    return prunedH

def rulesFromConseq(freqSet, H, supportData, br1, minConf=0.7):
    m = len(H[0])
    if (len(freqSet) > (m + 1)):
        Hmp1 = aprioriGen(H, m+1)
        Hmp1 = calcConf(freqSet, Hmp1, supportData, br1, minConf)
        if (len(Hmp1) > 1):
            rulesFromConseq(freqSet, Hmp1, supportData, br1, minConf)




if __name__ == "__main__":
   print "apriori 알고리즘"
   dataset = loadDataSet()
   L, suppData = apriori(dataset)
   print "L:"  + str(L)
   print "........................."
   print "suppData:" + str(suppData)

   rules = generateRules(L, suppData, minConf=0.7)

레퍼런스  : 머신러닝 인 액션 

요즘 데이터분석과 관련하여  텐서플로우와 스파크(ML)등의 머신러닝 솔루션들이 굉장히 유행하고 있습니다. 물론 저것들이 삶을 편안하게 만들어주기도 하지만 대부분의 데이터 분석은 저런 거창한 것 말고 평균,편차 같은 기본적인 개념으로 부터 시작되고 있으며 이러한 개념을 조금씩 변경해가며 더 의미있는 가치를 찾기 위해 빠르게 적용해보는 과정을 거쳐야하는데  그러기 위해서는 

1. 직접 코딩해서 기본적인 데이터분석 유틸리티 함수들을 만들어봐야한다. 

2. SQL문을 잘 다루어야한다. 

3. 엑셀을 잘 다루어야 한다. 

이 3가지는 기본이라고 저는 생각합니다. 소규모 데이터를 가지고 이리저리 반복해서 돌려보는 과정은 매우 중요하며 이런 기본적인 것들도 못하면서 하둡,텐서플로우나 깔짝대고, 데이터 분석 한다 라고 칭할 수는 없겠지요. (이것은 논쟁의 여지가 있습니다.) 그래서 이것들 중 1번에 대하여  "밑바닥부터 시작하는 데이터과학" 등의 좋은책의 내용을 통해서 살펴보는 시간을 갖도록 하겠습니다. 

통계,확률,패턴인식 분야의 내용들의 수식은 외우기가 힘듭니다. 외워도 금방 까먹어지는게 통계관련 공식인거 같습니다. (예를들어 정규분포,확률밀도라는 개념은 매우 쉽게 수긍이 가지만 , 그 식을 외우는건 좀 .;;) 또한 체득하고 나면 당연한거 아냐? 라고 느껴지는 알고리즘(수식) 인데 글로 설명하면 매우 산만해지는거 학문인거 같습니다. 하지만  우리들은 소프트웨어 개발자이기 때문에 수많은 기술변화도 따라가야하는 운명에 있는데 쓸때 없이 공식 및 수학을 외우고 있을 수는 없지 않겠습니까? 

제가 초등학교 5학년때 만(10,000) 단위 암산을 했었는데요. 주산,암산을 배워보신분은 아시겠지만 그것이 가능한것은 기억의 매개로 주판을 이용한 것 입니다. 마찬가지로 우리 개발자들은 그 기억의 매개체로 코드를 이용하면 개념과 함정(예를들어 K-Means 를 통해 군집화하면 길이차를 가지고 구분짓는 것이기 때문에 문제가 생기는 도메인 또한 많다) 에 대한 이해가 더 오래갈 것이라 보며 덤으로 가져다 사용도 할 수 있을것입니다.

* 그 상상의 매개체로써의 언어로 "파이썬" 을 선택했으며 정말 좋은 언어라고 생각합니다.

 * 참고로 "밑바닥부터 배우는 데이터과학" 서적은 numpy,scikit-learn 등의 외부라이브러리를 활용은 배제되었습니다.

K-Means 알고리즘 이란 ? 

머신러닝의 큰 바운더리로 중 하나로 군집을 들수있는데 K-Means 는 군집분석을 할때  활용되는 기본 알고리즘 입니다. 군집이란 , 쉽게 말해 도서관에서 여러분이 어떤 기준에 따라서 책을 그룹핑할때 하는 행동 . 그때 책들을 군집화 한다라고 말할수있습니다. (역사책, 수학책, 소프트웨어 책)  군집할때는 어떤 데이터간의 유사성을 정량적으로 측정하여 수치화하는게 핵심이며, (역사와 수학, 소프트웨어 를 각각 어떤 기준으로 수치화) 수치화가 된후 K-Means 알고리즘을 이용하여 자동 그룹핑 할 수 있게 되는데 K-Means 는 말 그래도 K 개로 그룹핑하라~ 와 같습니다. 유사한 입력 값끼리 묶어서 군집을 찾습니다.굉장히 클리어한 알고리즘이며 간단한데에 비해 많은곳에서 주요 알고리즘으로 사용됩니다.

다음 코드는 K-Means 를 나타내는데 그 수행과정을 보면 

1) 임의의 K개의 군집수를 결정하고 (여기선 3개) 그것의 초기치를 1개씩 할당(여기선 랜덤)하여 위치 설정한다.

self.means = random.sample(inputs, self.k)

2) 각각의 데이터에 대해 K개의 위치까지의 거리를 구하고 가장 가까운 군집에 소속시킨다.(유클리드 거리를 이용)

   def classify(self, input):
        return min(range(self.k),
                   key=lambda i: squared_distance(input, self.means[i]))

3) 3개의 군집으로 나뉘어진 데이터를 가지고 새로운 중앙의 위치를 설정한다.

   # 현재 클러스터링의 중심점 재 계산
    for i in range(self.k):
        i_points = [p for p, a in zip(inputs, assignments) if a == i] # i 번 요소들의 리스트
          if i_points:
            self.means[i] = vector_mean(i_points) # 그것들의 평균 ( 중심점 )

4) 이전 군집과 새롭게 중앙점을 이동한후에 얻어진 군집이 갖으면 종료한다.

        while True:
            # 요소들이 어느 (변경된)중심점에 가까이 있는지 매핑 [1,1,0,2,1, ... 이런식
            new_assignments = map(self.classify, inputs)

            # 기존 매핑과 동일하다면 더 이상 프로세싱 할 필요 없음. 리턴 ~
            if assignments == new_assignments:
                return

이 과정을 통하여 K개의 군집으로 나눈다. 아래는 전체 소스이다.

# coding: utf-8

from __future__ import division
from linear_algebra import squared_distance, vector_mean
import random


class KMeans:

    def __init__(self, k):
        self.k = k  # 클러스터의 갯수
        self.means = None  # 클러스터의 평균들

    # input 값이 3개의 means 중 어디와 가장 가까운지 찾아서 0,1,2 중 하나를 리턴
    def classify(self, input):
        return min(range(self.k),
                   key=lambda i: squared_distance(input, self.means[i]))

    def train(self, inputs):

        self.means = random.sample(inputs, self.k)
        assignments = None

        while True:
            # 요소들이 어느 (변경된)중심점에 가까이 있는지 매핑 [1,1,0,2,1, ... 이런식
            new_assignments = map(self.classify, inputs)

            # 기존 매핑과 동일하다면 더 이상 프로세싱 할 필요 없음. 리턴 ~
            if assignments == new_assignments:
                return

            # 현재 매핑을 저장
            assignments = new_assignments

            # 현재 클러스터링의 중심점 재 계산
            for i in range(self.k):
                i_points = [p for p, a in zip(inputs, assignments) if a == i] # i 번 요소들의 리스트
                if i_points:
                    self.means[i] = vector_mean(i_points) # 그것들의 평균 ( 중심점 )




if __name__ == "__main__":

    inputs = [[-14,-5],[13,13],[20,23],[-19,-11],[-9,-16],[21,27],[-49,15],[26,13],[-46,5],[-34,-1],[11,15],[-49,0],[-22,-16],[19,28],[-12,-8],[-13,-19],[-41,8],[-11,-6],[-25,-9],[-18,-3]]

    random.seed(0) # so you get the same results as me
    clusterer = KMeans(3)
    clusterer.train(inputs)
    print "3-means:"
    print clusterer.means
    print

문제 발생 1 

위처럼 단순 길이차이로 군집하는 알고리즘을 이용해서 

5,6,50,60,500,600 <-- 이걸 3개로 군집화 한다면 

[5,6,50,60] [500] [600]   이렇게 3개로 나뉠것이다.

하지만 [5,6] [50,60] [500,600]  이렇게 3개로 나뉘길 바랬더라면 망하는거다. OTL 

상향식 계층 군집화 

1. 각 데이터 포인트를 하나의 군집으로 간주한다.

2. 군집이 두개 이상이라면, 가장 가까운 두개의 군집을 찾아 하나의 군집으로 묶는다. 

위의 알고리즘으로 

5,6,50,60,500,600 <-- 이걸 3개로 군집화 한다면 

5에서 가장 가까운 6을 찾아서 [5,6] 군집화 

50에서 가장 가까운 60 찾아서 [50,60] 군집화

500 에서 가장 가까운 600 찾아서 [500,600] 군집화 

이렇게 나뉠것이다.

이런식으로 최소 거리를 이용한 계층 군집화는 위의  K-means 와는 다른 군집을 이룬다. 

K 값 찾기 

위 에서는 k 값을 그냥 3개로 했지만 k 값을 보다 합리적으로 찾을 수 있는 방법도 여러가지이다.
쉬운 방법은 k 값에 대해 중심점과 각 데이터 포인트 사이의 거리의 제곱합을 그래프로 그리고, 그 그래프가 어디서 꺽이는지 관찰하는 것이다.

def squared_clustering_errors(inputs, k):
    clusterer = KMeans(k)
    clusterer.train(inputs)
    means = clusterer.means
    assignments = map(clusterer.classify, inputs)

    return sum(squared_distance(input, means[cluster])
               for input, cluster in zip(inputs, assignments))

아래 처럼 꺽이는데 K=3 이 적절했던 거 같다.

전체소스

https://github.com/Insight-book/data-science-from-scratch/blob/master/code/ch19_clustering.py

분산환경에서 K-Means 알고리즘 알아보기 

[하마 인사이드] 1. HAMA with K-Means