세미나 참가자를 위한 Hidden Markov Model 예제

안녕하세요? 김태경 교수입니다. 오랫만에 글을 씁니다.

세미나 자료를 준비하다보니 학생들이 HMM을 알고 있나 하는 생각이 퍼뜩들어 자료를 적다가 다 지웠습니다... 대신에 그냥 바로 쓸 수 있도록 샘플이나 만들어보자는 생각에서 포스트를 올립니다.

우선...

Hidden Markov Model?

머신러닝 방법 중 하나입니다. 아마 요즘은 딥러닝에 밀려서 텍스트 예측 분야에서는 좀 안 쓰는 것 같습니다만... 사회과학 분야에서는 그렇지 않을 것 같습니다. 훗. 논문의 리뷰어 중 한 명도 "그냥 딥러닝 하지?" 이랬지만서도.

지진파나 음성, 주식 거래량, 필기체 분석 등에서 흔한 방법론입니다.

1960년대에 Baum 등이 은닉 계층을 마코브 모형에 추가했지요. 자세한 것은 Rabiner(1989) 논문을 참조하세요.

Markov Chain

확률론 시간에 배웠지요? 행렬의 성분이 확률값으로 이루어진 행렬입니다. 한 상태에서 다른 상태로 변할 확률이 먼 과거의 자취보다 바로 직전의 상태에 따라 결정된다는 것입니다. 흔히 상태-전이 행렬이라고도 합니다. 시장점유율이나 고객 상태 전이, 쿠폰의 사용이나 광고 회피(제가 지금 다루는 것) 등에서 실증 데이터를 바탕으로 한 상태-전이 행렬을 구성할 수 있습니다.

왜 Markov Chain?

아무 곳에서나 쓰라는 것은 아니구요... 시간적인 선후 관계가 있는 데이터를 다룰 때 필요합니다. 체중의 변화나 주식의 변화, 학습이나 글쓰기 등이 그렇지요. 운동 데이터도 마찬가지로 시간의 선후가 있습니다. 골기퍼가 바로 골을 넣어주기도 하지만 적어도 골대에서 골대로 가는 것은 아니니까요. 착실하게 시간을 들여 단계(혹은 상태)를 밟아야 합니다.

관측 벡터

순차 데이터는 관측(observation) 벡터로 표현합니다.

유한 시간 T 동안 관측된 데이터는 위와 같습니다. 관측의 전후는 서로 맞물려 있으니 다음 관측치가 나타날 확률은 조건부 확률로 표현될 것입니다. 이전 관측치 모두가 데이터로 사용된다고요? 아니죠. 마코브 체인이니 직전만 의미가 있습니다. 물론 2차, 3차 마코브 체인을 만들어서 그 전 시간 정보도 반영할 수는 있습니다. 모델이 복잡해져 문제죠.

관측치의 확률을 계산하는 것은 그냥 곱하기죠? 넘어가겠습니다.

은닉 상태-전이

만약 관찰된 것이 사실상 우리가 모르는 또 다른 상태 전이에 영향을 받는다면 어떨까요? 관찰만으로 은닉된 상태-전이의 확률분포를 알아낼 수 있을까요?

예를 들어 여러분이 어떤 감옥에 갇혀 있다고 합시다. 창문에는 빛 하나 들어오지 않아 바깥의 날씨를 알 수가 없습니다. 날씨를 아는 유일한 방법(관측)은 간수가 우산을 가지고 있는지 혹은 그렇지 않은지 뿐입니다. 그리고 용기를 내어 열 번 정도 간수에게 날씨가 어떤지를 물어봤다고 합시다. 그 대가는 혹독했지만 간신히 정보를 얻었다고 합시다.

P(우산=1|맑은 날)=0.1
P(우산=1|비)=0.8
P(우산=1|흐림)=0.3

더 이상 죄수는 간수에게 물어볼 용기가 안납니다(많이 맞았나보다). 이제 우산을 가져왔는지(=1), 안 가져 왔는지(=0) 만 가지고 날씨를 예측해야 합니다. 날씨 상태를 q라고 하고, 역시 관측을 o라고 하면, 관측되는 조건부 확률은 다음과 같습니다.

Naive Bayes와 마찬가지로 독립적인 관측열을 제낀다면, 우도는 다음과 같이:

만약, 1차 마코브 체이이면,

자 이렇게 계산하면 됩니다.

문제는.... 현실적으로 우도 계산은 다이내믹 프로그램으로 풀어야 하는데요... 실제 데이터로는 우도 계산하는 일이 장난이 아니게 됩니다. 학생들 골탕 먹이려면 Baum-Welch 알고리즘 알려주고 R로 풀어오라고 시키면 되겠습니다.

depmixS4로 풀기

R의 dependent mixture models를 풀기 위해 고안된 depmixS4로 골치아픈 fitting 문제를 해결하도록 하겠습니다. 아마 Netzer et al.(1997)의 마케팅 논문에서는 GAUSS를 썼던 것 같안데요... 후덜덜... 저는 depmixS4를 씁니다. 이번 사례는 Daniel의 Bob과 Alice 문제를 참고 했습니다.

Bob과 Alice의 카지노 문제입니다. 여러분은 고객 담당 매니저입니다. 밥과 앨리스는 카지노에 가서 돈을 많이 씁니다. 그런데 돈을 둘 다 100만원 단위로 쓴다고 합시다. 평균적으로 밥은 1200만원, 앨리스는 400만원을 씁니다. 여러분의 임무는 두 고객에게 각기 다른 서비스를 제공해야 한다는 것입니다. 그런데 문제는 이 둘은 큰 손들이라서 노출되기를 꺼려해서 사실 여러분은 누가 방에서 게임을 하는 지 알수가 없습니다. 주변에 물어봐도 알려주지도 않습니다. 오로지 카지노의 배팅 금액만을 알 수 있는데 이 정보로 현재 게임을 하는 사람이 밥인지 앨리스인지를 알아내야 합니다.

밥과 앨리스는 다음과 같은 종류의 게임을 합니다. 주사위 두개를 던져 합이 4 이상이 나오면 이깁니다. 배팅은 한 번만 할 수 있고 승자가 합니다. 만약 주사위 승부에서 진다면 상대방에게 찬스가 넘어갑니다.

라이브러리

필요한 라이브러리를 불러옵니다.

library(depmixS4) #HMM을 위해

library(ggplot2) #그래프를 그리려고

library(gridExtra) #그래프 배치를 위해

library(reshape2) #melt함수 쓰려고

library(dplyr) #연쇄 규칙 써먹으려고

시뮬레이션

우선 시뮬레이터를 만듭시다.

simulate<-function(N,dice.val=6,jbns,switch.val=4)

{

  bob.dice<-sample(1:dice.val,N,replace=T)+sample(1:dice.val,N,replace=T)

  alice.dice<-sample(1:dice.val,N,replace=T)+sample(1:dice.val,N,replace=T)

  bob.jbns<-rpois(N,jbns[1])

  alice.jbns<-rpois(N,jbns[2])

  #states

  draws<-data.frame(

    state=rep(NA,N),

    obs=rep(NA,N),

    dice=rep(NA,N))

  draws$state[1]<-'alice'

  draws$obs<-alice.jbns[1]

  draws$dice<-alice.dice[1]

  for(k in 2:N)

  {

    if(draws$state[k-1]=='alice')

    {

      if(draws$dice[k-1]<switch.val+1) {

        draws$state[k]<-'bob'

        draws$obs[k]<-bob.jbns[k]

        draws$dice[k]<-bob.dice[k]

      } else {

        draws$state[k]<-'alice'

        draws$obs[k]<-alice.jbns[k]

        draws$dice[k]<-alice.dice[k]

      }

    } else if(draws$state[k-1]=='bob') 

    {

      if(draws$dice[k-1]<switch.val+1) {

        draws$state[k]<-'alice'

        draws$obs[k]<-alice.jbns[k]

        draws$dice[k]<-alice.dice[k]

      } else {

        draws$state[k]<-'bob'

        draws$obs[k]<-bob.jbns[k]

        draws$dice[k]<-bob.dice[k]

      }

    }

  }

  return(cbind(roll=1:N,draws))

}

코드는 어렵지 않습니다. 그냥 밥과 앨리스의 승패를 가상으로 만들어 본 것입니다. rpois() 함수가 사용되었는데 상황에 따라 적절히 써야겠지요? 우리는 배팅 금액(100만원 단위)이니까 nonnegative 카운트 데이터입니다. 그래서 포아송 분포를 가정합니다.

이제 시뮬레이션을 수행합니다.

set.seed(201904007)

N=100

draws<-simulate(N,jbns=c(12,4),switch.val=4)

HMM Fitting

depmixS4가 있기 때문에 이후 일은 더 쉽습니다.

fit.hmm<-function(draws) {

  #HMM

  mod<-depmix(obs~1,

              data=draws,

              nstates=2,

              family=poisson())

  fit.mod<-fit(mod)

  #사후확률 계산합니다.

  est.states<-posterior(fit.mod)

  #결과를 정리하겠습니다.

  tbl<-table(est.states$state,draws$state) #상태도 확인

#앨리스? 혹은 밥? 

  whoswho<-c(colnames(tbl)[which.max(tbl[1,])],colnames(tbl)[which.max(tbl[2,])])

#머신러닝 결과값 추가합니다.

  draws$est.state.labels<-whoswho[est.states$state]

#정확도를 계산해봅시다. 

  confusion_matrix<-table(draws[,c('state','est.state.labels')])

  accuracy<-sum(diag(confusion_matrix))/sum(confusion_matrix)

#그래프를 그리기위해 데이터를 업데이트 합니다.

  est.states$roll<-1:100

  colnames(est.states)[2:3]<-whoswho

  hmm.post.df<-melt(est.states,measure.vars=c('alice','bob'))

#결과값을 돌려줍니다.

  return(list(draws=draws,hmm.post.df=hmm.post.df))

}

추정을 해 보겠습니다.

hmm1<-fit.hmm(draws)

그래프로 결과 확인하기

Daniel이 작성한 그래프 코드를 조금 고쳐서 쓰겠습니다. 4개의 그래프를 그리는데 위에서부터 차례대로, 1)실제 관측치, 2)실제 누구? 3)추정된 누구? 4)사후 확률값
이렇습니다. 느낌은 딥러닝에서 softmax 구하고 핫코딩 결과 보는 것 비슷합니다.

코드입니다.

mycols<-c("blue","green")
plot.hmm.output<-function(model.output) {
  g0<-(ggplot(model.output$draws,aes(x=roll,y=obs))+geom_line()+
         theme(axis.ticks=element_blank(),axis.title.y=element_blank()))%>%ggplotGrob()
  g1<-(ggplot(model.output$draws,aes(x=roll,y=state,fill=state,col=state))+
         geom_bar(stat='identity',alpha=I(0.7))+
         scale_fill_manual(values=mycols,
                           name='State:\nPerson\nthat\nrolled\ndice',
                           labels=c('Alice','Bob'))+
         scale_color_manual(values=mycols,
                            name='State:\nPerson\nthat\nrolled\ndice',
                            labels=c('Alice','Bob'))+
         theme(axis.ticks=element_blank(),
               axis.text.y=element_blank())+
         labs(y='Actual state'))%>%ggplotGrob()
  g2<-(ggplot(model.output$draws,aes(x=roll,y=est.state.labels,
                                     fill=est.state.labels,
                                     col=est.state.labels))+
         geom_bar(stat='identity',alpha=I(0.7))+
         scale_fill_manual(values=mycols,
                           name='State:\nPerson that\nrolled the\ndice',
                           labels=c('Alice','Bob'))+
         scale_color_manual(values=mycols,
                            name='State:\nPerson that\nrolled the\ndice',
                            labels=c('Alice','Bob'))+
         theme(axis.ticks=element_blank(),
               axis.text.y=element_blank())+
         labs(y='Estimated State'))%>%ggplotGrob()
  g3<-(ggplot(model.output$hmm.post.df,
              aes(x=roll,y=value,col=variable))+
         geom_line()+
         scale_color_manual(values=mycols,
                            name='State:\nPerson that\nrolled the\ndice',
                            labels=c('Alice','Bob'))+
         theme(axis.ticks=element_blank(),
               axis.text.y=element_blank())+
         labs(y='Posterior Prob.'))%>%ggplotGrob()
  g0$widths<-g1$widths
  return(grid.arrange(g0,g1,g2,g3,widths=1,nrow=4))
}

이제 결과를 보겠습니다.

plot.hmm.output(hmm1)

그래프는 다음과 같습니다.

마무리

HMM은 관측되는 확률 특성과 은닉된 확률 특성으로 모델링할 수 있습니다. Netzer 등은 이러한 인사이트로 상당히 괜찮은 아이디어를 선보였습니다. 직접적인 효과(direct effect)로 인한 관측된 확률 특성, 그리고 간접적 효과(indirect effect)로 인한 은닉된 확률 특성을 바탕으로 어떤 정책이나 디자인의 효과에 따른 차이를 시뮬레이션 할 수 있다는  것입니다. 이때 직접적 효과 특성에 특이성을 반영함으로써 개인 수준, 그룹 수준, 조직 수준 등으로 범위를 바꿀 수 있습니다.

시간에 따라 변하는 은닉 마코브 체인 특성은 정보시스템 연구에서 얼마든지 찾을 수 있습니다. 게임에서 레벨의 변화, 혹은 MMORPG에서 역할(마법사, 탱커, 암살자, 원딜, 서폿) 변화, 고객의 등급(잠재군, 진입군, 안정군, 위험군, 이탈군)이나 광고에 대한 태도(부정적, 중립적, 긍정적), 혹은 자기 통제 수준(낮은 수준, 중간 수준, 높은 수준) 등 다양한 간접적 이유들이 존재합니다. 이와 같은 특성들에 따라 각종 전략적 행동을 취할 수 있다는 점도 생각해야겠습니다.

추가하자면...

논문 제출 후 어떤 MIS Quarterly 리뷰어의 요청에 따라 (ㅜㅜㅜ) 오토인코딩+RNN으로 접근해 본 적도 있는데요... 데이터가 무지막지하게 필요했습니다...만 데이터를 제공하는 쪽에서 완전 난색을 표현해서 무산된 바가 있습니다. 아하... 회사 입장도 이해가 되고... 여러 노력을 펼쳤지만 App 데이터를 전면적으로 받기 위해 사용자 동의를 일일이 받는 것도 부담이 되는지라...

쉬운 일이 없네요.