Rによる隠れマルコフモデル

大学院のゼミでPRMLを読んでいる。私の担当は13章「系列データ」の後半の「状態空間モデル」（SSM: State Space Model)である。13章の前半は「隠れマルコフモデル」(HMM: Hidden Markov Model)について記述されている。SSMは日本のお家芸であり、統計数理研究所の先代所長の北川先生や現所長の樋口先生により発展した。ところが、SSMの離散バージョンだと言われているHMMはなぜかあまり日本ではクローズアップされていない。PRMLはHMM→SSMの順番で記述されているので、丁度良い機会だからついでにHMMについて勉強しようと思ってこのブログをアップした。


とりあえず、はじめからアルゴリズムを実装するのは私の実力では厳しいので、Rのパッケージで遊んでみよう思ってあれこれ探したところ、{RHmm}というパッケージが見つかった。このパッケージについてはいくつかの日本語ブログで言及されているものの、既にCRANから削除されてしまっている。なんじゃ。と思い、いろいろ調べてみると、{depmixS4}というHMMのパッケージが見つかった。これは最近でも更新されているし、R-bloggerでも紹介されているので間違いないと思って、こちらをいじってみることにした。
（本記事では理論については書いていない。理論について知りたい方はPRML等を参照のこと。）


このパッケージにはSpeedというサンプルデータがあって、これを基にサンプルコードが展開されている。ヘルプには詳しく書かれていなかったが、どうやら心理学のスピードと正確性のトレードオフに関する実験データのようである。


■Speedデータ

rt: レスポンス時間(response time)
corr: 正答スコア(1/0) (accuracy scores)
pacc: 正確性に注意を払っている度合い (the pay-off for accuracy)
prev: 1期前の正答スコア(1/0) (accuracy scores (0/1) on the previous trial)

■Rコード
次のコードによりパッケージを読み込んでSpeedデータを可視化する。

#####Library
library(depmixS4)
#####Data
data(speed)
#####Visualization
plot(as.ts(speed[1:168,]), main = "Speed-accuracy trade-off")

次に隠れマルコフモデルを構築する。手始めに隠れ状態 = 2 と想定する。隠れ状態は恣意的に設定するが、AICやBICを参照することで最適な隠れ状態を設定することが可能であるっぽい（勉強中）。以下のコードにて、depmixオブジェクトを作って、EMアルゴリズムでパラメータを推定。

##### Modeling
set.seed(1)
mod <- depmix(response = rt ~ 1, data = speed, nstates = 2, trstart = runif(4))
fm  <- fit(mod) #EMアルゴリズムによるパラメータ推定

次にモデル推定結果を出力。先ずは情報量の確認。

print(fm)

# Convergence info: Log likelihood converged to within tol. (relative change) 
# 'log Lik.' -88.73058 (df=7)
# AIC:  191.4612 
# BIC:  220.0527 

ちなみに、隠れ状態 = 1 の場合には、AIC = 614.6619, BIC = 622.8309 となることから、上記の隠れ状態 = 2 の方がAICもBICも低いため良いモデルという評価となった。次に推定されたパラメータを確認。

summary(fm)

# Initial state probabilties model 
# pr1 pr2 
# 1   0 
# 
# 
# Transition matrix 
# toS1       toS2
# fromS1 0.9156683 0.08433171
# fromS2 0.1164891 0.88351091
# 
# Response parameters 
# Resp 1 : gaussian 
# Re1.(Intercept)    Re1.sd
# St1        6.385053 0.2441601
# St2        5.510363 0.1917513

Transition matrixを見ると、推移確率は以下のように推定されたということ。
S1⇒S1: 0.9156683
S1⇒S2: 0.08433171
S2⇒S2: 0.1164891
S2⇒S1: 0.88351091

次に、"Response parameters"の解釈だが、これはおそらく推定された各隠れ状態(S1・S2）におけるレスポンス時間(rt)の平均と分散を記述していると思われる（勉強中のため断言できず）。

最後に、Viterbiアルゴリズムによって潜在変数系列を推定した上で結果を可視化すると次のようになる。

##### Viterbiアルゴリズムによる系列推定
result <- posterior(fm)

result$rt <- speed$rt
result$time <- seq(1,439)

#rt-state
gg1 <- ggplot(data=result,aes(x=time, y=rt)) +
        geom_rect(aes(fill = state, xmin = time - 0.5, xmax = time + 0.5, ymin=4.5,
                    ymax = max(rt) + 0.1)) +
        geom_line(colour="red", size=0.5)
print(gg1)

#corr-state
gg2 <- ggplot(data=result,aes(x=time, y=corr)) +
        geom_rect(aes(fill = state, xmin = time - 0.5, xmax = time + 0.5, ymin=1.0,
                      ymax = 2.0)) +
        geom_line(colour="red", size=0.5)
print(gg2)

レスポンス時間が落ちていて、かつ誤答が多い部分が認識されていることが分かります。このように隠れマルコフモデルを用いることで、時系列データを離散分離することができるようです。未だ勉強始めたばかりでわからないことだらけですが、そのうちまた続編を書きたいと思います。

■参考文献

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

• 作者: C.M. ビショップ
• 出版社/メーカー: 丸善出版
• 発売日: 2012/02/29
• メディア: 単行本

やっぱりこれが日本語で読める一番詳しい本です。状態空間モデルとの関係が分かりやすいです。

Hidden Markov Models for Time Series: An Introduction Using R (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

• 作者: Walter Zucchini
• 出版社/メーカー: Chapman and Hall/CRC
• 発売日: 2009/04/30
• メディア: ハードカバー

洋書ですがRコードが開示されていて分かりやすいです。具体的事例もいくつかあります。

タグ：隠れマルコフモデル Hidden Markov model R HMM depmixS4

CentOS6.5にRStudio Serverをインストール

WindowsXPのサポート終了に伴い、自宅PCのOSをCentOSに移行しました。なので、分析環境もいちから出直し。Rのインストールはすんなりいったのだけど、RStudio Serverがインストールできない。頑張って解決したので、Tipsをメモっておきます。

■インストール

rpm -Uvh rstudio-server-0.98.507-i686.rpm

ところが、これではインストールできない。

エラー: 依存性の欠如:
libR.so は rstudio-server-0.98.507-1.i686 に必要とされています
libRblas.so は rstudio-server-0.98.507-1.i686 に必要とされています
libRlapack.so は rstudio-server-0.98.507-1.i686 に必要とされています

→Errorになってしまった。失敗。

■やり直し

rpm -Uvh --nodeps rstudio-server-0.98.507-i686.rpm

→こっちで成功。"--nodeps"というのがポイントだったみたい。

■以降のコマンド
・参照サイト
http://d.hatena.ne.jp/hiratake55/20120527/1338091311

・インストールできているか確認

rstudio-server verify-installation

・RStudio Serverを起動

rstudio-server start

・システム起動時に自動起動する

chkconfig rstudio-server on
chkconfig --list rstudio-server

DVD付 CentOS徹底入門 第3版

• 作者: 中島 能和
• 出版社/メーカー: 翔泳社
• 発売日: 2012/03/09
• メディア: 大型本

Rで学ぶデータ・プログラミング入門 ―RStudioを活用する―

• 作者: 石田 基広
• 出版社/メーカー: 共立出版
• 発売日: 2012/10/24
• メディア: 単行本

タグ：RStudio R CentOS linux

ROC曲線をggplot2で描画

ROC曲線をRの美麗な可視化パッケージであるggplot2で表示するプログラムを作ってみました。

ROCを描画するパッケージとしてEpiというものがあるのだが、ちょっと出力される画像が好みではないのでggplot2でやってみました。三重大学の奥村先生のサイトを参考にしました。

ROC <- function(score, actual) {
  ###ROC-AUC計算
  o = order(score, decreasing = T) 
  fp = tp = fp_prev = tp_prev = 0
  nF = sum(actual == F)
  nT = sum(actual == T)
  score_prev = -Inf
  ber_min = Inf
  area = 0
  rx = ry = numeric(length(o))
  n = 0
  
  #cat("i", "\t", "n", "\t", "score[j]", "\t", "tp", "\t", "fp", "\t", "rx[n]", "\t", "ry[n]", "\t", "area", "\t", "AUC", "\n")
  
  for (i in seq_along(o)) {
    j = o[i]
    if (score[j] != score_prev) {
      area = area + (fp - fp_prev) * (tp + tp_prev) / 2
      n = n + 1
      rx[n] = fp/nF
      ry[n] = tp/nT
      ber = (fp/nF + 1 - tp/nT)/2
      AUC = area/(nF*nT) 
      
      #cat(i, "\t", n, "\t", score[j], "\t", tp, "\t", fp, "\t", rx[n], "\t", ry[n], "\t", area, "\t", AUC, "\n")
      
      if (ber < ber_min) {
        ber_min = ber
        th = score_prev
        rx_best = fp/nF
        ry_best = tp/nT
      }
      score_prev = score[j]
      fp_prev = fp
      tp_prev = tp
    }
    if (actual[j] == T) {
      tp = tp + 1
    } else {
      fp = fp + 1
    }
  }
  area = area + (fp - fp_prev) * (tp + tp_prev) / 2
  AUC = area/(nF*nT)
  n = n + 1
  rx[n] = fp/nF  # = 1
  ry[n] = tp/nT  # = 1
  
  ###ggplot2で可視化
  AUC <- round(AUC, digits = 4)
  tmp <- data.frame(x=rx, y=ry)
  label <- paste("AUC = ", AUC)
  gg <- ggplot(tmp) + geom_path(aes(x,y, color = "red"), size = 1.0) + 
    geom_path(data=data.frame(x = c(0,1), y = c(0,1)),aes(x,y), colour = "gray", size = 1.0) +
    scale_x_continuous("False_Positive Rate (1 - Specificity)", labels=percent, limits = c(0, 1)) +
    scale_y_continuous("True_Positive Rate (Sensivity or Recall)", labels=percent,limits = c(0, 1)) +
    theme(legend.position = "none") +
    ggtitle("ROC Curve") +
    geom_text(data = NULL, x = 0.85, y = 0.15, label = label, colour = "black")
  print(gg)
  
  ###ROC-AUC計算
  cat("\n")
  cat("AUC =", area/(nF*nT), "th =", th, "\n")
  cat("BER =", (rx_best + (1-ry_best))/2,
      "OR =", (ry_best/(1-ry_best))/(rx_best/(1-rx_best)), "\n")
  print(table(score >= th, actual, dnn=c("Predicted","Actual")))
}

####実行例
library(ggplot2)
library(scales) #ggplot2のscaleを%表示にするのに必要

s = c(16,15,14,13,12,11,10, 9, 8, 8, 8, 8, 7, 6, 5) #score
a = c( T, T, F, T, T, T, F, T, T, T, T, F, F, T, F) #actual
ROC(s, a)

【実行結果】

やっぱりggplot2はキレイでいいですね！

◆参考文献
・マシンラーニング (Rで学ぶデータサイエンス 6)

マシンラーニング (Rで学ぶデータサイエンス 6)

• 作者: 辻谷 将明
• 出版社/メーカー: 共立出版
• 発売日: 2009/06/09
• メディア: 単行本

ROCやAUCの基本的な考え方はこの本の１章にあります。

タグ：R ggplot2 Roc AUC データマイニング