Data Science by R and Python

統計学を、広く、深く、わかりやすく。

Hidden Markov Model - by R

隠れマルコフモデルをつくってみる。

こんにちは、日曜日も終わりにさしかかっておりますが、今日は完全に息抜きをしたくて、HMMでも勉強して、Rで書いてみるかということで、作成してみました。

f:id:tomoshige_n:20150510180834p:plain

潜在変数がスイッチになって、観測値を発生させる分布が変わるというものなので、Rのコードとしては、こんな感じです。

library(MASS)

transition = matrix(c(9/10,0,4/10,1/10,7/10,0,0,3/10,6/10),3,3)
switch = function(z,transition){
	trans_prob = transition[z,]
	next_state = sample(c(1,2,3),1,prob=trans_prob)
	return(next_state)
}

sample_nmix = function(z){
	r = 0
	if(z == 1){
		VCM = matrix(c(0.5,0.25,0.5,0.25),2,2)
		r = mvrnorm(1,mu=c(2,3),Sigma=VCM)
	}
	if(z == 2){
		VCM = matrix(c(0.5,0.25,0.5,0.25),2,2)
		r = mvrnorm(1,mu=c(-3,0),Sigma=VCM)
	}
	if(z == 3){
		VCM = matrix(c(0.5,0.25,0.5,0.25),2,2)
		r = mvrnorm(1,mu=c(2,-3),Sigma=VCM)
	}
	return(r)
}

iter = 100
class.r = rep(0,iter)
sample.r =data.frame(x = rep(0,iter), y = rep(0,iter))
z_initial = which(rmultinom(1,1,prob=c(1/3,1/3,1/3))==1)
x_initial = sample_nmix(z_initial)
class.r[1] = z_initial
sample.r[1,] = x_initial

for(i in 2:iter){
	class = switch(class.r[i-1],transition)
	class.r[i] = class
	sample.r[i,] = sample_nmix(class)
}


plot(sample.r,col=1,type="l",lwd=0.5)
points(sample.r,col=class.r,type="p")

Rodeo!

そういえば、昨日Yhatが出している、Rodeoを少しだけ触ってみましたが、良い感じです!github.com

→今回のもRodeoで作れば良かった...

参考文献

言わずと知れた...

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

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