【求助统计牛人】如何用R做HMM模型选择？

小妹在纠结毕业论文中……
建立了几个隐马尔可夫模型，想选出一个最优的，据说可以基于BIC准则之类的做模型选择

问题是：R的HMM包里没有带模型选择功能

不知R里有没有做模型选择的通用程序（无论基于什么准则都行），或者有高人可以帮忙编个程么……

给个思路也行啊……bow～～

Matlab会将iteration过程记录下来
以上述例子为例,
iteration=13时,达到converged.
logliks也由 -554.4527逐渐增加到 -454.1974(maximum likelihood )
这时 logL 就取 -454.1974

parameters的算法
请参考范例page 1-3
em_notes.pdf
em_notes.pdf (82.81 KB)

2011-5-15 16:21:56 上传

number of states : N=3
number of symbols: M=4={e,f,g,h}
parameters = N+(N-1)*N+(N-1)*M =17

Hidden Markov Models--Model selection and checking
ch6_slides.pdf (280.28 KB)

2011-5-12 11:49:11 上传

建议楼主使用package"HiddenMarkov"
function 较齐全
function logLik() 可算出 Log Likelihood of Hidden Markov Model
然后依公式就可算出 AIC or BIC

1. Akaike's information criteria: AIC = -2log L + 2p.
2. Bayesian information criteria: BIC = -2log L + plogN,
     where L is the likelihood of the fitted model,
     p is the number of parameters,
     and N is the number of data points.

1. >install.packages("HiddenMarkov")
   
2. >require("HiddenMarkov")
   
3. >help("HiddenMarkov")

复制代码

试一试...

2# epoh
楼上能不能将全部的讲义share出来啊？

2# epoh
谢谢指点。但HiddenMarkov包里貌似需要假定分布的？想想也是哦，求最大似然当然需要明确是什么分布了……

问题是，我是用一串序列建立的HMM，不大好假定分布的，这样一来还能有办法进行模型选择么…？

ch6_slides系引用
Hidden Markov Models for Time Series An Introduction Using R
这本书的教材.
大家有兴趣的话,可上作者网站
下载R code,data,及各章节slides

package "HMM" 采用
observation,transition probabilities,and emission probabilities
方法和 matlab相同.
然matlab有function  hmmtrain()
[guessTR,guessE,logliks] = hmmtrain(seqs,guessTR,guessE,varargin)
可以算出 log likelihood
如果楼主有用matlab的话,可以试试
在R的话,需要用到
observation,bw$hmm$transProbs,bw$hmm$emissionProbs加以编程

######in R
library(HMM)
# Initial HMM
hmm = initHMM(c("A","B"),c("L","R"),
transProbs=matrix(c(.9,.1,.1,.9),2),
emissionProbs=matrix(c(.5,.51,.5,.49),2))
print(hmm)
# Sequence of observation
a = sample(c(rep("L",100),rep("R",300)))
b = sample(c(rep("L",300),rep("R",100)))
observation = c(a,b)
# Baum-Welch
bw = baumWelch(hmm,observation,10)
print(bw$hmm)
#$transProbs
#    to
#from            A           B
#   A 9.974870e-01 0.002513041
#   B 5.296986e-06 0.999994703
#$emissionProbs
#      symbols
#states         L         R
#     A 0.2485981 0.7514019
#     B 0.7497916 0.2502084


%%%%%%%in Matlab
tr=[0.9 0.1; 0.1 0.9]
e=[0.50 0.50; 0.51 0.49]
load('123_HMM.txt')
x=X123_HMM';
[guessTR,guessE,logliks] = hmmtrain(x,tr,e)  

guessTR =

   0.9975   0.0025
   0.0000   1


guessE =

   0.2486   0.7514
   0.7498   0.2502

logliks=-554.4527 -553.5753 -546.9715 -518.1480 -484.5900
              -471.5483 -463.1453 -457.1819 -454.6034 -454.2221
              -454.1987 -454.1975 -454.1974

6# epoh
多谢高人指点，我安装了Matlab，也跑了你给我的程序，但是看不懂结果……为什么Loglikes会出来这么多个结果呢？AIC和BIC准则里的L怎么取呢……？以及对不假定分布的HMM来说，准则里的p值取多少呢？求教……

2# epoh
不好意思喔，因为是第一次发悬赏帖，我不知道选出最佳答案后100个币是不是会自动到你账上？如果你没收到，就随便上传个附件设置成100个币出售吧，我来下载，谢谢～～ o(∩_∩)o

哈哈!楼主太客气了!
倒是这个问题,我无法回答.
因为我也不知道.
楼主不用在意!