読者です 読者をやめる 読者になる 読者になる

StatsFragments

Python, R, Rust, 統計, 機械学習とか

R の {MSwM} パッケージでマルコフ転換モデルをためす

マルコフ転換モデルとは

数式を使わない説明。

サーモンとインターネット広告とマルコフ転換モデル|インターネット広告代理店で働くデータサイエンティストのブログ

マルコフ状態転換モデルのRパッケージ{MSwM}の使い方(異常値検出・ステータス変化検出などに有用) - 銀座で働くData Scientistのブログ

Exampleとしては vignettes がわかりやすい。

http://cran.r-project.org/web/packages/MSwM/vignettes/examples.pdf

まず

手元には下図のような時系列データがあるとする。経験的になんとなく 目的変数 y は説明変数 x に関係していそう、ということがわかっている。が、ぱっと見でも y のレベルが途中でシフトしており、常に一定の関係がなりたっているわけではなさそう (背後に複数の状態がありそう)。

こんなときにマルコフ転換モデルを使うと 背後にある状態変化を踏まえて モデルを作ってくれる。

f:id:sinhrks:20141026195725p:plain

上のサンプルデータは以下のようにして作った。時点に応じて x, y の関係が切り替わるデータとなっている。

※説明変数は目的変数とまったく関係のないものでもよいが、(一部の状態にでも)関係していそうなものを選んだほうがよい。説明変数が目的変数と関係ないものだと 切片の変化 = ベースラインの変化でしか状態転換を検出できないと思う。その場合マルコフ転換を使う必要がない。

x <- rpois(500, lambda = 10)   # x はポアソン分布に従う
y1 <- x * 4 + 20    # 状態1のときの x, yの関係
y2 <- x * 2 + 60    # 状態2のときの x, yの関係

# ノイズ付与
noise <- rnorm(1:500, mean = 10, sd = 5)
y1 <- y1 + noise
y2 <- y2 + noise

# 適当なタイミングで状態転換 (1:200までは状態1, 201:400までは状態2, 以降は状態1)
y <- c(y1[1:200], y2[201:400], y1[401:500])

# 実際に観測されたデータ
observed <- data.frame(x = x, y = y)

x と y1, y2の関係を図示すると下図のようになる。x, yに二つの状態があることがわかっていれば、 x, yはこのように見える。

状態別に分離したデータ

f:id:sinhrks:20141026201256p:plain

が、観測者にとっては 状態がそもそもいくつあるのかすらわかっていない。そのため x, yは 状態の区別なくひとまとまりで見える。マルコフ転換モデルでは、観測データがいくつかの状態から発生していると捉えて、うまいこと上図のように分離してくれる。

観測されるデータ

f:id:sinhrks:20141026215845p:plain

マルコフ転換モデルの作成

{MSwM} でマルコフ転換モデルを作るのは非常にカンタンで、

  1. 目的変数/説明変数の関係を 線形回帰モデル, もしくは一般化線形モデルとして記述する
  2. 1で作ったモデル式とパラメータを MSwM::msmFit に渡してマルコフ転換モデルをフィットさせる

だけ。

モデル式作成

ここでは y と x は線形関係にあるものとして、lm 関数でモデル式を作る。この状態 (全体をひとつの状態として線形回帰した場合)で summary 表示して Coefficients をみると、ちょうど y1, y2の間くらいに線を引いていることがわかる。残差が大きいので 調整済みR二乗値も低い。

library(MSwM)

m.lm <- lm(y ~ x, data = observed)
summary(m.lm)
# 
# Call:
# lm(formula = y ~ x, data = observed)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -24.303  -9.354  -1.914   9.617  29.224 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  45.7468     1.7202   26.59   <2e-16 ***
# x             3.2262     0.1636   19.71   <2e-16 ***
# ---
# Signif. codes:  
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 11.51 on 498 degrees of freedom
# Multiple R-squared:  0.4383, Adjusted R-squared:  0.4372 
# F-statistic: 388.7 on 1 and 498 DF,  p-value: < 2.2e-16

モデルのフィット

上で作ったモデルを msmFit に渡せばよい。パラメータの意味は、

  • k: マルコフ転換モデルの状態数。 ここでは 背後に2つの状態があるものとして指定。
  • sw: 各パラメータが状態転換時に変化するかどうかを Logical で指定
  • p: ARモデルの係数
  • family: (GLMの場合) 確率分布の family
m.mswm <- msmFit(m.lm, k = 2, sw = c(TRUE, TRUE, TRUE, TRUE), p = 1)
summary(m.mswm)
# Markov Switching Model
# 
# Call: msmFit(object = m.lm, k = 2, sw = c(TRUE, TRUE, TRUE, TRUE), 
#     p = 1)
# 
#        AIC      BIC    logLik
#   3038.846 3101.397 -1513.423
# 
# Coefficients:
# 
# Regime 1 
# ---------
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)(S)  69.3263     4.0606 17.0729   <2e-16 ***
# x(S)             2.1795     0.1187 18.3614   <2e-16 ***
# y_1(S)          -0.0103     0.0429 -0.2401   0.8103    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 4.99756
# Multiple R-squared: 0.6288
# 
# Standardized Residuals:
#           Min            Q1           Med            Q3           Max 
# -1.431396e+01 -2.056292e-02 -1.536781e-03 -1.098923e-05  1.584478e+01 
# 
# Regime 2 
# ---------
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)(S)  30.2820     1.7687 17.1210   <2e-16 ***
# x(S)             3.9964     0.0913 43.7722   <2e-16 ***
# y_1(S)          -0.0045     0.0203 -0.2217   0.8245    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 4.836684
# Multiple R-squared: 0.8663
# 
# Standardized Residuals:
#           Min            Q1           Med            Q3           Max 
# -13.202056966  -0.771854514   0.002211602   1.162769110  12.417873232 
# 
# Transition probabilities:
#             Regime 1    Regime 2
# Regime 1 0.994973376 0.003347279
# Regime 2 0.005026624 0.996652721

出力中の Regime 1, Regime 2が 背後にある(と考えられる)二つの状態をあらわす。y の真の値と比較すると、出力中の Regime 1が y2 <- x * 2 + 60 に対応しており、

# Regime 1 
# ---------
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)(S)  69.3263     4.0606 17.0729   <2e-16 ***
# x(S)             2.1795     0.1187 18.3614   <2e-16 ***
# y_1(S)          -0.0103     0.0429 -0.2401   0.8103   

Regime 2が y1 <- x * 4 + 20 と対応していることがわかる。

切片のずれが少し気になるがノイズのせいかな、、?調整済みR二乗値から全体として改善しているとは言える。

# Regime 2 
# ---------
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)(S)  30.2820     1.7687 17.1210   <2e-16 ***
# x(S)             3.9964     0.0913 43.7722   <2e-16 ***
# y_1(S)          -0.0045     0.0203 -0.2217   0.8245    

モデルのプロット

状態 (Regime) ごとに 目的変数 + 指定した説明変数と、その状態にある確率を陰影でプロット

plotReg(m.mswm, c('x'), regime = 1)

f:id:sinhrks:20141026210629p:plain

各時点がそれぞれどちらの状態 (Regime)にあるかを確率でプロット

plotProb(m.mswm, which = 1)

f:id:sinhrks:20141026210603p:plain

各時点の状態と変化点を取得

ある時点でどちらの状態 (Regime) にいるかを知りたい場合は、その時点で一番 発生している(していた)確率の高い Regimeを拾ってくればよい。

probable <- apply(m.mswm@Fit@smoProb[-1, ], 1, which.max)
> probable
  [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [30] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
...

異常値 / 変化点は Regime が転換した時点になるので、上で取得したベクトルの値が変化した箇所が TRUE になるようにすると、

c(FALSE, diff(probable) != 0)
  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [11] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
...
[181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[191] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
[201] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
...
[381] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[391] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
[401] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
...
[491] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

ということで 最初のデータ作成時に指定した変化点 (201, 401番目)に近いところが検知できているのがわかる。