Expectation–maximization과 R을 이용한 구현

/*

 * http://sosal.kr/
 * made by so_Sal
 */

Expectation Maximization 알고리즘에 대해 어렴풋이 알고있었는데,

이번기회에 R프로그래밍을 이용하여 직접 구현해보았다.

라이브러리를 사용하면 쉽게 사용할 수 있지만,

수리통계와 확률론에 평소에 약했다고 느끼는 터라, 이번 기회에 한번 직접 구현해보았다.

EM 알고리즘을 적용할 수 있는 예는 방대하지만, Clustering을 예제로 하여 이 글을 작성한다.

Clustering의 대상은 DNA Sequence로 한다.

샘플: 데이터

변수: 하나의 데이터(샘플)가 가지는 다양한 값

Clustering은 비슷한 변수를 가지고 있는 샘플들을 군집화 시켜주는 알고리즘이다.

K-means

Clustering 기법중, k-means clustering은 각 샘플이 가진 변수의 값들이

임의로 정한 군집의 중심과 얼마나 가깝냐에 따라서 샘플을 그룹화 하는 것을 기본으로 한다.

군집의 중심값을 그룹의 중심으로 계속 이동시키면서 샘플을 clustering 하는 기법이다.

군집의 중심값이 완벽히 Convergence 할 수 있고, 데이터가 너무 큰 경우 적당히 convergence 했을 때 멈출 수 있다.

기댓값 최대화 알고리즘(expectation-maximization algorithm)

Expectation-maximization 알고리즘은 모델 파라미터 (Model parameter)가 바로 하나의 군집이 된다.

각 데이터에 대해서, 모든 모델파라미터의 likelihood를 구할 수 있는데, 이때 최대 likelihood 값을 가지는 모델로 군집화 한다.

각 모델파라미터에 속하는 데이터들을 다시 모아서, 모델의 파라미터를 수정하는 방식으로 진행된다.

likelihood (가능도, 우도)

하나의 모델로부터 주어진 값이 나올 확률을 가능도라고 한다.

예를 들어, 첫 번째 동전에서 앞면이 나올 확률을 40%, 뒷면이 나올 확률을 60%라고 하고

두 번째 동전에서 앞면이 나올 확률을 60%, 뒷면이 나올 확률을 40%라고 하자.

동전을 던졌더니 '앞뒤앞앞' 가 나왔다. 

첫 번째 동전과 두번쨰 동전을 기반으로 한 likelihood를 구하면 다음과 같다.

첫 번째 동전의 likelihood: 40% * 60% * 40% * 40% = 3.84%

두 번째 동전의 likelihood: 60% * 40% * 60% * 60% = 8.64%

Model parameter

모델 파라미터는 위에서 언급한 첫번째 동전에서, 앞면이 나올 확률, 뒷면이 나올 확률이 된다.

따라서 '앞뒤앞앞' 이라는 샘플은 두 번째 동전의 모델 파라미터에서 만들어졌을 확률이 높다.

바로 이 likelihood가 k-means에서 중심값과의 distance라고 생각하면 된다.

k-means에서, 하나의 샘플로부터 가까운 군집을 선택하는것이, 

바로 EM에서는 하나의 샘플로부터 더 높은 확률의 likelihood에 해당하는 모델(군집)을 선택하는 것이 된다.

모든 동전 Sequence에서 이렇게 likelihood를 계산하여 어떤 모델로부터 만들어진 Sequence인지 선택을 한 후,

같은 모델로부터 만들어졌다고 평가받은 모든 Sequence를 기반으로 다시 모델 파라미터를 선정하는 것이 바로 EM 알고리즘이다.

Nature biotechnology에서 이 Expectation maximization 알고리즘을 쉽게 설명한 논문이 있다.

Do, Chuong B., and Serafim Batzoglou. "What is the expectation maximization algorithm?." Nature biotechnology 26.8 (2008): 897-899.

Figure 1은 아래와 같다.

Figure 1. Parameter estimation for complete and incomplete data. (a) Maximum likelihood estimation. For each set of ten tosses, the maximum likelihood procedure accumulates the counts of heads and tails for coins A and B separately. These counts are then used to estimate the coin biases. (b) Expectation maximization. 1. EM starts with an initial guess of the parameters. 2. In the E-step, a probability distribution over possible completions is computed using the current parameters. The counts shown in the table are the expected numbers of heads and tails according to this distribution. 3. In the M-step, new parameters are determined using the current completions. 4. After several repetitions of the E-step and M-step, the algorithm converges.

이 논문에서도 동전의 앞, 뒷면이 나온 Sequence를 이용하여 EM 알고리즘을 통해 어떤 동전으로부터 나온 Sequence인지 예측하는 예가 된다. 단순히 Sequence를 Clustering 하는 문제로 해석할 수 있으면서, 데이터로부터 모델링 후에 새로 들어온 sequence 데이터에 대해서 어떤 모델로부터 추론된 sequence인지 prediction 하는 모델로도 해석할 수 있다.

Figure 1a의 경우, 주어진 Sequence로부터 어떤 Model (Coin A, Coin B)에서 추론되었는지를 likelihood를 이용하여 계산하는 모습이다.

EM 알고리즘은 바로 이 Model (Coin A, Coin B)의 파라미터를 최적으로 하는 값을 찾고자 하는게 목표이다.

모델의 수는 2개 (Coin A, B) 이기 때문에, 달리 해석하면 주어진 5개의 동전(HT) Sequence를 가장 잘 clustering할 수 있는 최적의 2가지 모델 파라미터를 셋팅하는 것이다.

Figure 2b에 이 EM 알고리즘이 잘 나와있는데,

E-step과 M-step의 반복이 바로 EM 알고리즘의 원리이다.

1. Model parameter setting

주로 k-means에서 각 군집의 중심값을 무작위로 선정하는 것 처럼, EM 알고리즘에서도 각 모델 파라미터의 값을 무작위로 선정한다. 이 때, Gaussian distribution을 사용하는 것이 일반적이다.

2. E-step (Expectation)

설정된 각 모델들의 파라미터로부터, 각 샘플의 likelihood를 계산한다.

k개의 모델이 있다면, 하나의 샘플로부터 k개의 likelihood가 계산될 것이다.

하나의 샘플로부터 추론된 k개의 Likelihood로부터, 최대값이 되는 likelihood가 어떤 모델 파라미터로부터 만들어진 것인지 확인한다. 그 샘플은 최대값이 되는 모델로부터 추론되었을 것이라고 가정한다. (즉 그 샘플은 최대 likelihood를 구성하는 모델에 속한다고 가정한다.)

3. M-step (Maximization)

E-step을 진행했다는 뜻은, 각 모델에 속하는 샘플들이 어떻게 구성되어있는지 알고있다는 뜻이다.

따라서 각 모델에 속하는 sequence를 기반으로 Model의 parameter를 다시 계산한다.

4. Repeat!

likelihood값 혹은 parameter값이 바뀌지 않을 때 까지, 2번과 3번을 반복한다.

http://cse-wiki.unl.edu/wiki/index.php/Expectation_maximization

- R programming을 이용한 구현

리스트를 이용하여 구현하였다.

정확하게 구현한건지는.. 누가 검증해준 바가 없기 때문에..

(혹시 어색한 부분이 있다면 조언 부탁드립니다.)

Clustering 잘 되는지 결과를 보니 어느정도 올바르게 작동하는 것 같긴 하다.

가독성을 최대한 높이기 위해 노력을 했긴 했지만.. 쉽진 않았다.

처음에 예상했던것 보다는 코드가 짧게 구현된 것 같다 (?).

Objective function을 바꾸고자 한다면 아래 Main code에서 Likelihood 대신에 원하는 objective function을 적용하면 된다.

# Expectation function

Expectation <- function(Sequences, Families, nSequence, lSequence, nFamily) {

    pSequences <- list()

    for (i in 1:nSequence) {

        pSequences[[i]] <- matrix(0, length(Sequences[[i]]), nFamily)

        colnames(pSequences[[i]]) <- paste0("Family", 1:nFamily)

        rownames(pSequences[[i]]) <- 1:lSequence

        for (k in 1:nFamily)

            pSequences[[i]][, k] <- Families[[k]][cbind(1:lSequence, Sequences[[i]])]

        }

    return(pSequences)

}

 

# Maximization function

Maximization <- function(Families, Sequences, Cluster, lSequence, nFamily, Characters) {

    for (i in 1:nFamily) {

        FamilySequence <- Sequences[which(Cluster == i)]

        for (l in 1:lSequence) {

            pseudo <- c(unlist(lapply(FamilySequence, '[[', l)), Characters) # Characters: PseudoCount

            FamilySequenceArgument <- table(pseudo) # Model parameter가 0이 되는걸 방지

            Families[[i]][l,] <- FamilySequenceArgument / sum(FamilySequenceArgument)

        }

    }

    return(Families)

}

 

#0. Data setting.

lSequence <- 8

nSequence <- 10

Characters <- c("A", "C", "G", "T")

Sequences <- c("AGTGTAAC",

                "AGTTTATC",

                "CGATTTTC",

                "AGATTTAC",

                "AGATTATC",

                "GGGGGGCA",

                "CGCGGTGA",

                "GGGGGACT",

                "GGCGGGGT",

                "GGGTGGGT")

Sequences <- strsplit(Sequences, "")

 

#1. Hidden variable setting.

nFamily = 2

Families <- list()

for (i in 1:nFamily) {

    Families[[i]] <- matrix(runif(lSequence * length(Characters)), nrow = lSequence)

    Families[[i]] <- Families[[i]] / rowSums(Families[[i]])

    colnames(Families[[i]]) <- Characters

    rownames(Families[[i]]) <- 1:lSequence

}

######### Main #########

maxIter = 100

for (n in 1:maxIter) {

    #2. E-step

    pSequences <- Expectation(Sequences, Families, nSequence, lSequence, nFamily)

    pSequences_log <- lapply(pSequences, log)

    logLikelihood <- lapply(pSequences_log, colSums)

    Cluster <- c()

    for (i in 1:length(logLikelihood))

        Cluster[i] <- which(logLikelihood[[i]] == max(logLikelihood[[i]]))

    # Objective function: Log-likelihood의 합.

    # Likelihood를 곱하는것 보다, log를 취한 값을 더하는 것이 더 효율적.

 

    #3. M-step.

    Families <- Maximization(Families, Sequences, Cluster, lSequence, nFamily, Characters)

}