HOWTOUSETHEHMMTOOLBOX(MATLAB)
一、离散输出的隐马尔科夫模型(DHMM,HMMwithdiscreteoutputs)
最大似然参数估计EM(BaumWelch算法)
Thescriptdhmm_em_anexampleofhowtolearnanHMMwith
rebeQ=2statesandO=terandom
stochasticmatricesasfollows.
O=3;
Q=2;
prior0=normali(rand(Q,1));
transmat0=mk_stochastic(rand(Q,Q));
obsmat0=mk_stochastic(rand(Q,O));
Nowwesamplenex=20quencesoflengthT=10eachfromthismodel,touas
trainingdata.
T=10;%序列长度
nex=20;%样本序列数目
data=dhmm_sample(prior0,transmat0,obsmat0,nex,T);
akearandomguessastowhattheparametersare,
prior1=normali(rand(Q,1));%初始状态概率
transmat1=mk_stochastic(rand(Q,Q));%初始状态转移矩阵
obsmat1=mk_stochastic(rand(Q,O));%初始观测状态到隐藏状态间的概率转移矩
阵
andimproveourguessusing5iterationsofEM...
[LL,prior2,transmat2,obsmat2]=dhmm_em(data,prior1,transmat1,obsmat1,
'max_iter',5);
%prior2,transmat2,obsmat2为训练好后的初始概率,状态转移矩阵及混合状态
概率转移矩阵
LL(t)isthelog-likelihoodafteriterationt,sowecanplotthelearningcurve.
序列分类
Toevaluatethelog-likelihoodofatrainedmodelgiventestdata,proceedasfollows:
loglik=dhmm_logprob(data,prior2,transmat2,obsmat2)%HMM测试
Note:thediscretealphabetisassumedtobe{1,2,...,O},whereO=size(obsmat,2).
Hencedatacannotcontainany0s.
Toclassifyaquenceintooneofkclass,trainupkHMMs,oneperclass,andthen
computethelog-likelihoodthateachmodelgivestothetestquence;ifthei'th
modelisthemostlikely,thendeclaretheclassofthequencetobeclassi.
Computingthemostprobablequence(Viterbi)
FirstyouneedtoevaluateB(i,t)=P(y_t|Q_t=i)forallt,i:
B=multinomial_prob(data,obsmat);
Thenyoucanu
[path]=viterbi_path(prior,transmat,B)
二、具有高斯混合输出的隐马尔科夫模型(GHMM,HMMwithmixtureof
Gaussiansoutputs)
MaximumlikelihoodparameterestimationusingEM(BaumWelch)
Letusgeneratenex=50vector-valuedquencesoflengthT=50;eachvectorhassize
O=2.
O=2;
T=50;
nex=50;
data=randn(O,T,nex);
NowletufitamixtureofM=2GaussiansforeachoftheQ=2statesusing
K-means.
M=2;
Q=2;
left_right=0;
prior0=normali(rand(Q,1));
transmat0=mk_stochastic(rand(Q,Q));
[mu0,Sigma0]=mixgauss_init(Q*M,reshape(data,[OT*nex]),cov_type);
mu0=reshape(mu0,[OQM]);
Sigma0=reshape(Sigma0,[OOQM]);
mixmat0=mk_stochastic(rand(Q,M));
Finally,letusimprovetheparameterestimatesusingEM.
[LL,prior1,transmat1,mu1,Sigma1,mixmat1]=...
mhmm_em(data,prior0,transmat0,mu0,Sigma0,mixmat0,'max_iter',2);
SinceEMonlyfindsalocaloptimum,tialisation
procedureillustratedaboveisverycrude,andisprobablynotadequateforreal
applications...Clickhereforareal-worldexampleofEMwithmixturesofGaussians
usingBNT.
Whattodoifthelog-likelihoodbecomespositive?
Itispossibleforp(x)>1ifp(x)isaprobabilitydensityfunction,suchasa
Gaussian.(Therequirementsforadensityarep(x)>0forallxandint_xp(x)=1.)In
practicethisusuallymeansyourcovarianceisshrinkingtoapoint/deltafunction,so
youshouldincreathewidthoftheprior(ebelow),orconstrainthematrixtobe
sphericalordiagonal,orclampittoalargefixedconstant(notlearnitatall).Itisalso
veryhelpfultoensurethecomponentsofthedatavectorshavesmallandcomparable
magnitudes(ue.g.,KPMstats/standardize).
Thisisawell-knownpathologyofmaximumlikelihoodestimationforGaussian
mixtures:theglobaloptimummayplaceonemixturecomponentonasingledata
point,andgiveit0covariance,allyreliesonthe
factthatEMcannotfindtheglobaloptimumtoavoidsuchpathologies.
Whattodoifthelog-likelihooddecreasduringEM?
SinceIimplicitlyaddapriortoeverycovariancematrix(ebelow),whatincreas
isloglik+log(prior),butwhatIprintisjustloglik,whichmayoccasionallydecrea.
T
betterinitializationorfewerclusters(states).
Whattodoifthecovariancematrixbecomessingular?
Estimatesofthecovariancematrixoftenbecomesingularifyouhavetoolittledata,
oriftoofewpointsareassignedtoaclustercenterduetoabadinitializationofthe
ca,youshouldconstrainthecovariancetobesphericalordiagonal,
oradjusttheprior(ebelow),ortryabetterinitialization.
HowdoIaddapriortothecovariancematrix?
BuriedinsideofKPMstats/mixgauss_Mstepyouwillethatcov_priorisinitialized
to0.01*
changethis,youwillneedtomodifythemhmm_emfunctionsoitcalls
mixgauss_Mstepwithadifferentvalue.
Sequenceclassification
Toclassifyaquence(e.g.,ofspeech)intooneofkclass(e.g.,thedigits0-9),
proceedasintheDHMMcaabove,bututhefollowingproceduretocompute
likelihood:
loglik=mhmm_logprob(data,prior,transmat,mu,Sigma,mixmat);
Computingthemostprobablequence(Viterbi)
FirstyouneedtoevaluateB(t,i)=P(y_t|Q_t=i)forallt,i:
B=mixgauss_prob(data(:,:,ex),mu,Sigma,mixmat);
wheredata(:,:,ex)y,u
[path]=viterbi_path(prior,transmat,B);
三、具有高斯输出的HMM
ThisisjustlikethemixtureofGaussiansca,exceptwehaveM=1,andhencethere
isnomixingmatrix.
OnlineEMfordiscreteHMMs/POMDPs
Forsomeapplications(e.g.,reinforcementlearning/adaptivecontrol),itisnecessary
iptdhmm_em_online_demogivesanexampleofhow
todothis.
隐马尔科夫(HMM)实例分析
A和B是好朋友,但是他们离得比较远,每天都是通过电话了解对方那天作了什么.B
仅仅对三种活动感兴趣:公园散步,购物以及清理房间.他选择做什么事情只凭当天天气.A对
于B所住的地方的天气情况并不了解,但是知道总的趋势.在B告诉A每天所做的事情基础
上,A想要猜测B所在地的天气情况.
A认为天气的运行就像一个马尔可夫链.其有两个状态“雨”和”晴”,但是无法直接观察
它们,也就是说,它们对于A是隐藏的.每天,B有一定的概率进行下列活动:”散步”,“购物”,或
“清理”.因为B会告诉A他的活动,所以这些活动就是A的观察数据.这整个系统就是一个隐
马尔可夫模型HMM.
A知道这个地区的总的天气趋势,并且平时知道B会做的事情.也就是说这个隐马尔可夫
模型的参数是已知的.下面是概率转移矩阵和两种天气下各种活动的概率:
雨天晴天
雨天0.70.3
晴天0.40.6
散步购物清理
雨天0.10.40.5
晴天0.60.30.1
下面是一段程序来描述各个变量。
//状态数目,两个状态:雨或晴
states=(‘Rainy’,‘Sunny’)
//每个状态下可能的观察值
obrvations=(‘walk’,‘shop’,‘clean’)
//初始状态空间的概率分布
start_probability={‘Rainy’:0.6,‘Sunny’:0.4}
//与时间无关的状态转移概率矩阵
transition_probability={
’Rainy’:{‘Rainy’:0.7,‘Sunny’:0.3},
’Sunny’:{‘Rainy’:0.4,‘Sunny’:0.6},
}
//给定状态下,观察值概率分布,发射概率
emission_probability={
’Rainy’:{‘walk’:0.1,‘shop’:0.4,‘clean’:0.5},
’Sunny’:{‘walk’:0.6,‘shop’:0.3,‘clean’:0.1},
}
在这些代码中,start_probability代表了A对于B第一次给她打电话时的天气情况的不
确定性(A知道的只是那个地方平均起来下雨多些).在这里,这个特定的概率分布并非平衡的,
平衡概率应该接近(在给定变迁概率的情况下){‘Rainy’:0.571,‘Sunny’:0.429}。
transition_probability表示马尔可夫链下的天气变迁情况,在这个例子中,如果今天下雨,
那么明天天晴的概率只有30%.代码emission_probability表示了B每天作某件事的概率.
如果下雨,有50%的概率他在清理房间;如果天晴,则有60%的概率他在外头散步。
A和B通了三天电话后发现第一天B去散步了,第二天他去购物了,第三天他清理房
间了。A现在有两个问题:这个观察序列“散步、购物、清理”的总的概率是多少?(注:这
个问题对应于HMM的基本问题之一:已知HMM模型λ及观察序列O,如何计算P(O|λ)?)
最能解释这个观察序列的状态序列(晴/雨)又是什么?(注:这个问题对应HMM基本问
题之二:给定观察序列O=O1,O2,…OT以及模型λ,如何选择一个对应的状态序列S=
q1,q2,…qT,使得S能够最为合理的解释观察序列O?)
至于HMM的基本问题之三:如何调整模型参数,使得P(O|λ)最大?这个问题事实上就
是给出很多个观察序列值,来训练以上几个参数的问题。
**************************
最大似然参数估计EM(BaumWelch算法)
Thescriptdhmm_em_anexampleofhowtolearnanHMMwithdiscreteoutputs.
LettherebeQ=2statesandO=terandomstochasticmatricesas
follows.
O=3;
Q=2;
prior0=normali(rand(Q,1));
transmat0=mk_stochastic(rand(Q,Q));
obsmat0=mk_stochastic(rand(Q,O));
Nowwesamplenex=20quencesoflengthT=10eachfromthismodel,touastrainingdata.
T=10;%序列长度
nex=20;%样本序列数目
data=dhmm_sample(prior0,transmat0,obsmat0,nex,T);
akearandomguessastowhattheparametersare,
prior1=normali(rand(Q,1));%初始状态概率
transmat1=mk_stochastic(rand(Q,Q));%初始状态转移矩阵
obsmat1=mk_stochastic(rand(Q,O));%初始观测状态到隐藏状态间的概率转移矩阵
andimproveourguessusing5iterationsofEM...
[LL,prior2,transmat2,obsmat2]=dhmm_em(data,prior1,transmat1,obsmat1,'max_iter',5);
%prior2,transmat2,obsmat2为训练好后的初始概率,状态转移矩阵及混合状态概率转移矩
阵
LL(t)isthelog-likelihoodafteriterationt,sowecanplotthelearningcurve.
序列分类
Toevaluatethelog-likelihoodofatrainedmodelgiventestdata,proceedasfollows:
loglik=dhmm_logprob(data,prior2,transmat2,obsmat2)%HMM测试
Note:thediscretealphabetisassumedtobe{1,2,...,O},whereO=size(obsmat,2).Hencedata
cannotcontainany0s.
Toclassifyaquenceintooneofkclass,trainupkHMMs,oneperclass,andthencompute
thelog-likelihoodthateachmodelgivestothetestquence;ifthei'thmodelisthemostlikely,
thendeclaretheclassofthequencetobeclassi.
Computingthemostprobablequence(Viterbi)
FirstyouneedtoevaluateB(i,t)=P(y_t|Q_t=i)forallt,i:
B=multinomial_prob(data,obsmat);
Thenyoucanu
[path]=viterbi_path(prior,transmat,B)
二、具有高斯混合输出的隐马尔科夫模型(GHMM,HMMwithmixtureofGaussians
outputs)
MaximumlikelihoodparameterestimationusingEM(BaumWelch)
Letusgeneratenex=50vector-valuedquencesoflengthT=50;eachvectorhassizeO=2.
O=2;
T=50;
nex=50;
data=randn(O,T,nex);
NowletufitamixtureofM=2GaussiansforeachoftheQ=2statesusingK-means.
M=2;
Q=2;
left_right=0;
prior0=normali(rand(Q,1));
transmat0=mk_stochastic(rand(Q,Q));
[mu0,Sigma0]=mixgauss_init(Q*M,reshape(data,[OT*nex]),cov_type);
mu0=reshape(mu0,[OQM]);
Sigma0=reshape(Sigma0,[OOQM]);
mixmat0=mk_stochastic(rand(Q,M));
Finally,letusimprovetheparameterestimatesusingEM.
[LL,prior1,transmat1,mu1,Sigma1,mixmat1]=...
mhmm_em(data,prior0,transmat0,mu0,Sigma0,mixmat0,'max_iter',2);
SinceEMonlyfindsalocaloptimum,tialisationprocedure
illustratedaboveisverycrude,andisprobablynotadequateforreal
applications...Clickhereforareal-worldexampleofEMwithmixturesofGaussiansusingBNT.
Whattodoifthelog-likelihoodbecomespositive?
Itispossibleforp(x)>1ifp(x)isaprobabilitydensityfunction,suchasaGaussian.(The
requirementsforadensityarep(x)>0forallxandint_xp(x)=1.)Inpracticethisusuallymeans
yourcovarianceisshrinkingtoapoint/deltafunction,soyoushouldincreathewidthoftheprior
(ebelow),orconstrainthematrixtobesphericalordiagonal,orclampittoalargefixed
constant(notlearnitatall).Itisalsoveryhelpfultoensurethecomponentsofthedatavectors
havesmallandcomparablemagnitudes(ue.g.,KPMstats/standardize).
Thisisawell-knownpathologyofmaximumlikelihoodestimationforGaussianmixtures:the
globaloptimummayplaceonemixturecomponentonasingledatapoint,andgiveit0covariance,
allyreliesonthefactthatEMcannotfindtheglobal
optimumtoavoidsuchpathologies.
Whattodoifthelog-likelihooddecreasduringEM?
SinceIimplicitlyaddapriortoeverycovariancematrix(ebelow),whatincreasisloglik+
log(prior),butwhatIprintisjustloglik,ggeststhatone
tterinitializationorfewerclusters
(states).
Whattodoifthecovariancematrixbecomessingular?
Estimatesofthecovariancematrixoftenbecomesingularifyouhavetoolittledata,oriftoofew
poica,you
shouldconstrainthecovariancetobesphericalordiagonal,oradjusttheprior(ebelow),ortrya
betterinitialization.
HowdoIaddapriortothecovariancematrix?
BuriedinsideofKPMstats/mixgauss_Mstepyouwillethatcov_priorisinitializedto0.01*I.
gethis,youwill
needtomodifythemhmm_emfunctionsoitcallsmixgauss_Mstepwithadifferentvalue.
Sequenceclassification
Toclassifyaquence(e.g.,ofspeech)intooneofkclass(e.g.,thedigits0-9),proceedasinthe
DHMMcaabove,bututhefollowingproceduretocomputelikelihood:
loglik=mhmm_logprob(data,prior,transmat,mu,Sigma,mixmat);
Computingthemostprobablequence(Viterbi)
FirstyouneedtoevaluateB(t,i)=P(y_t|Q_t=i)forallt,i:
B=mixgauss_prob(data(:,:,ex),mu,Sigma,mixmat);
wheredata(:,:,ex)y,u
[path]=viterbi_path(prior,transmat,B);
三、具有高斯输出的HMM
ThisisjustlikethemixtureofGaussiansca,exceptwehaveM=1,andhencethereisnomixing
matrix.
OnlineEMfordiscreteHMMs/POMDPs
Forsomeapplications(e.g.,reinforcementlearning/adaptivecontrol),itisnecessarytolearna
iptdhmm_em_online_demogivesanexampleofhowtodothis.
HMM有三个典型(canonical)问题:
1、已知模型参数,计算某一特定输出序列的概率.通常使用forward算法解决.
2、已知模型参数,寻找最可能的能产生某一特定输出序列的隐含状态的序列.通常使用
Viterbi算法解决.
3、已知输出序列,寻找最可能的状态转移以及输出概率.通常使用Baum-Welch算法以及
Viterbi算法解决.Reverd
另外,最近的一些方法使用Junctiontree算法来解决这三个问题。
x=importdata('output_')%读出output_此文件中的数据
变量说明:
设有M个状态,N个符号Markov模型。
TRANS:对应状态转移矩阵,大小为M*M,表示各状态相互转换的概率,
TRANS(i,j)表示从状态i转换到状态j的概率。
EMIS:对应符号生成矩阵,又叫混淆矩阵,
观察符号概率分布。EMIS(i,j)代表在状态i时,产生符号j的概率。
函数介绍:
hmmgenerate—GeneratesaquenceofstatesandemissionsfromaMarkovmodel
从一个马尔科夫模型产生状态序列和输出序列,该序列具有模型所表达的随机性特征。
Arandomquenceqofemissionsymbols
Arandomquencestatesofstates
用法:
[q,states]=hmmgenerate(len,TRANS,EMIS)
hmmgenerate(...,'Symbols',SYMBOLS)
hmmgenerate(...,'Statenames',STATENAMES)
示例:
trans=[0.95,0.05;0.10,0.90];%M=2
emis=[1/61/61/61/61/61/6;1/101/101/101/101/101/2];%N=6
[q,states]=hmmgenerate(100,trans,emis)%q=1*100,1,2,3,4,5,6其中的一个,因
为%eims=2*6,状态states为两个,所以1*100,1,2其中之一,因为状态转移矩阵为2*2
[q,states]=hmmgenerate(100,trans,emis,'Symbols',{'one','two','three','four','five','six'},...
'Statenames',{'fair';'loaded'})
%q,states分别用不同的字符来表示,将数字改成了字符
估计状态序列:EstimatingtheStateSequence
hmmviterbi—CalculatesthemostprobablestatepathforahiddenMarkovmodel
GiventhetransitionandemissionmatricesTRANSandEMIS,thefunctionhmmviterbiusthe
Viterbialgorithmtocomputethemostlikelyquenceofstatesthemodelwouldgothroughto
generateagivenquenceqofemissions:
给定状态转移矩阵TRANS和混淆矩阵EMIS,hmmviterbi使用Viterbi算法计算该模型最相似
的状态序列。
用法:
STATES=hmmviterbi(q,TRANS,EMIS)
hmmviterbi(...,'Symbols',SYMBOLS)
hmmviterbi(...,'Statenames',STATENAMES)
示例:
likelystates=hmmviterbi(q,TRANS,EMIS);
trans=[0.95,0.05;
0.10,0.90];
emis=[1/61/61/61/61/61/6;
1/101/101/101/101/101/2];
[q,states]=hmmgenerate(100,trans,emis);
estimatedStates=hmmviterbi(q,trans,emis);%由100个观察序列以及HMM矩阵(已知状
态转移矩阵和混淆矩阵,即确定了一个HMM模型)来求状态序列,此时状态用数字(1,2)
来表示
[q,states]=hmmgenerate(100,trans,emis,'Statenames',{'fair';'loaded'});
estimatesStates=hmmviterbi(q,trans,emis,'Statenames',{'fair';'loaded'});
%由100个观察序列以及HMM矩阵(已知状态转移矩阵和混淆矩阵,即确定了一个HMM
模型)来求状态序列,此时状态用状态名字来表示(fairorloaded)来表示
hmmestimate和hmmtrain用于通过给定的输出序列估计转移矩阵TRANS和混淆矩阵EMIS。
hmmestimate—Calculatesmaximumlikelihoodestimatesoftransitionandemission
probabilitiesfromaquenceofemissionsandaknownquenceofstates
通过一个输出序列和已知的状态序列,计算转移概率和输出概率的最大似然估计。
用法:
[TRANS,EMIS]=hmmestimate(q,states)
hmmestimate(...,'Symbols',SYMBOLS)
hmmestimate(...,'Statenames',STATENAMES)
hmmestimate(...,'Pudoemissions',PSEUDOE)
hmmestimate(...,'Pudotransitions',PSEUDOTR)
示例:
通过已知的输出序列和状态序列估计出转移状态和混淆矩阵。
[TRANS_EST,EMIS_EST]=hmmestimate(q,states)
TRANS_EST=
0.89890.1011
0.05850.9415
EMIS_EST=
0.17210.17210.17490.16120.18030.1393
0.58360.07410.08040.07890.07260.1104
hmmtrain—Calculatesmaximumlikelihoodestimatesoftransitionandemissionprobabilities
fromaquenceofemissions.在知道输出序列,不知道状态转移序列,但是对转移矩阵和混
淆矩阵有个初始猜测的情况下,可以使用hmmtrain估计转移状态和混淆矩阵,改函数可以
选择使用BaumWelch或者Viterbi算法,通过迭代的方式进行估计,可以设置迭代次数
Maxiterations和精度Tolerance。
用法:
[ESTTR,ESTEMIT]=hmmtrain(q,TRGUESS,EMITGUESS)
hmmtrain(...,'Algorithm',algorithm)
hmmtrain(...,'Symbols',SYMBOLS)
hmmtrain(...,'Tolerance',tol)
hmmtrain(...,'Maxiterations',maxiter)
hmmtrain(...,'Verbo',true)
hmmtrain(...,'Pudoemissions',PSEUDOE)
hmmtrain(...,'Pudotransitions',PSEUDOTR)
示例:
TRANS_GUESS=[.85.15;.1.9];
EMIS_GUESS=[.17.16.17.16.17.17;.6.08.08.08.0808];
YouestimateTRANSandEMISasfollows:
[TRANS_EST2,EMIS_EST2]=hmmtrain(q,TRANS_GUESS,EMIS_GUESS)
TRANS_EST2=
0.22860.7714
0.00320.9968
EMIS_EST2=
0.14360.23480.18370.19630.23500.0066
0.43550.10890.11440.10820.11090.1220
hmmdecode—Calculatestheposteriorstateprobabilitiesofaquenceofemissions
用于估计后验状态概率。对于一个输出序列q的后验状态概率,是指这个模型在确定的状
态下产生q中一个符号的条件概率。
PSTATES=hmmdecode(q,TRANS,EMIS)
TheoutputPSTATESisanM-by-Lmatrix,whereMisthenumberofstatesandListhelengthof
S(i,j)istheconditionalprobabilitythatthemodelisinstateiwhenitgeneratesthejth
symbolofq,giventhatqimitted.
PSTATES是一个M×L的矩阵,M是状态的个数,L是输出序列q的长度。PSTATES(i,j)是在
状态i时,产生q的第j个符号的条件概率。
细心的读者可能会发现,上面都没有使用到Markov模型的初始状态概率矩阵。实际上,如
果指定初始状态的概率,则上述函数默认从第一个状态开始,即初始状态为第一个状态的概
率是1。但是statisticstoolbox也提供了修改初始状态概率矩阵的方法。
设实际的状态转移矩阵和混淆矩阵分别是TRANS和EMIS,初始状态的概率分布为
p=[p1,p2,p3,...,pm],则可以通过以下方式重新设置初始状态。
TRANS_HAT=[0p;zeros(size(TRANS,1),1)TRANS];
EMIS_HAT=[zeros(1,size(EMIS,2));EMIS];
隐马尔科夫模型HMM
发表于2013年9月18日由了凡春秋USTC
目录
马尔科夫链
隐马尔科夫模型(HMM)
HMM简介
分析隐马尔科夫模型
产生测试序列
估计状态序列
估计状态转移矩阵和输出矩阵
估计后验状态概率
改变初始状态分布
马尔科夫链
如下图的马尔科夫模型状态转移矩阵
图中的方框代表你要建模过程的可能状态,箭头表示状态间的转移情况。箭头上的标签
表示转移概率。在过程中的每一步,模型会产生一个输出(决定于所处的状态)然后转移到
另一个状态。马尔科夫模型的一个重要特性是下一个状态只与当前的状态有关,与过去的转
移无关。
例如,对于一个顺序抛掷硬币的过程,存在两个状态:正和反。最近一次的投掷决定了
模型当前的状态,下一次投掷决定了如何转移到下一状态。如果硬币是公平的,则转移概率
都是二分之一,输出就是当前的状态。在更复杂的模型中,每一个状态可以以一个随机过程
产生输出,如可以通过投掷骰子来决定输出。
马尔科夫链是离散状态马尔科夫模型的数学描述。马尔科夫链如下定义:
一系列状态{1,2,…,M};
一个MXM的状态转移矩阵T,其中Tij表示从状态i到j的状态转移概率。矩阵每行的
和为1,因为它表示从一个给定状态到其他各状态的所有转移概率的和。
一系列可能的输出{s1,s2,…,sN},默认时可以当做{1,2,…,N},N为可能输出的个数,当
然你可以选择不同的数据集或符号来表示;
一个MXN的输出矩阵E,其中Eik表示模型在状态i时的输出sk的概率。
马尔科夫链的过程:从第0步的初始状态i0开始;然后以概率T1i1转移到状态i1,
并以概率Ei1k1输出sk1;最后,以概率T1i1Ei1k1Ti1i2Ei2k2…Tir-1irEirk
观察到状态i1i2…ir,同时输出sk1sk2…skr。
隐马尔科夫模型(HMM)
HMM简介
隐马尔科夫模型是你只观察到输出,但不知道模型产生输出所经历的状态。隐马尔科夫
模型分析就是从观察的数据恢复相关的状态序列。
例如,考虑一个有2个状态6个可能输出的马尔科夫模型,模型中用到
一个红骰子,6个面,分别标号1~6;
一个绿骰子,12个面,5个面标号2~6,其他面标记1;
一个非均衡的红硬币,以0.9的概率掷出正,0.1的概率掷出反;
一个非均衡的绿硬币,以0.95的概率掷出正,0.05的概率掷出反。
模型以下面的规则产生一系列来自{1,2,3,4,5,6}的数字序列:
先掷红骰子,记下骰子的输出;
然后掷红硬币:
如果是正,则掷红骰子记下结果;
否则掷绿骰子记结果;
接下来的下一步,掷与上一步所用骰子同颜色的硬币,如果是正,则掷同颜色的骰子;
否则掷另一个骰子。
模型的状态转移图如下
这个模型中,通过掷与状态同颜色的骰子来决定输出,通过掷与状态同颜色的硬币来决
定转移到哪一状态。
状态转移矩阵为
输出矩阵为
这个模型是非隐的,因为你从硬币、骰子的颜色可以知道状态序列。如果假设有人产生
了一系列输出,但不给你展示硬币和骰子,所有你能看到的只是输出序列。如果你开始看到
1出现的次数比其他都多,你可能会假设模型处于绿状态,但你不能确定,因为你不能看到
所掷骰子的颜色。
隐马尔科夫模型提出了如下问题:
1)给定输出序列,求最可能的状态序列?
2)给定输出序列,你如何估计模型的转移概率和输出概率?
3)模型得出给定序列的先验概率?
4)模型在序列某一位置是特定状态的后验概率?Hmmdecode
分析隐马尔科夫模型
与隐马尔科夫模型有关的统计工具箱函数如下
Hmmgenerate——从一个马尔科夫模型产生状态、输出序列;
Hmmestimate——从输出序列和已知的状态序列计算转移、输出序列概率的极大似然估
计;
Hmmtrain——从输出序列计算转移、输出序列概率的极大似然估计;
Hmmviterbi——计算隐马尔科夫模型的最大可能状态序列;
Hmmdecode——计算输出序列的状态后验概率。
下面来讲如何使用这些函数分析隐马尔科夫模型。
产生测试序列
下面的命令创建状态转移矩阵和输出矩阵
TRANS=[.9.1;.05.95;];
EMIS=[1/6,1/6,1/6,1/6,1/6,1/6;...
7/12,1/12,1/12,1/12,1/12,1/12];
要产生一组状态和输出的随机序列,使用函数hmmgenerate
[q,states]=hmmgenerate(1000,TRANS,EMIS);
输出参数q为输出序列,states为状态序列。
Hmmgenerate函数第0步开始于状态1,然后第一步转移到状态i1,而i1作为states
的第一个元素。如果要改变初始状态,看后面的改变初始状态分布一节。
估计状态序列
给定状态转移矩阵TRANS和输出矩阵EMIS,函数hmmviterbi使用Viterbi算法计算产
生给定输出序列q的极大可能状态序列
likelystates=hmmviterbi(q,TRANS,EMIS);
likelystates与q长度相同。
为了测试hmmviterbi的精确度,计算实际状态序列states与估计状态序列
likelystates吻合的百分比
sum(states==likelystates)/1000
ans=
0.8200
在这个例子中,估计的状态和实际吻合达82%
估计状态转移矩阵和输出矩阵
函数hmmestimate和hmtrain可以用来根据给定的输出序列估计状态转移矩阵和输出矩
阵。
使用hmmestimate
此函数需要你知道产生输出序列的状态序列。基于输出序列和状态序列,如下估计转移
矩阵和输出矩阵
[TRANS_EST,EMIS_EST]=hmmestimate(q,states)
TRANS_EST=
0.89890.1011
0.05850.9415
EMIS_EST=
0.17210.17210.17490.16120.18030.1393
0.58360.07410.08040.07890.07260.1104
与原始的序列对比
TRANS
TRANS=
0.90000.1000
0.05000.9500
EMIS
EMIS=
0.16670.16670.16670.16670.16670.1667
0.58330.08330.08330.08330.08330.0833
使用hmmtrain
如果你不知道状态序列,但知道TRANS和EMIS矩阵的初始猜测,可以使用hmmtrain
函数估计。
假设已知TRANS和EMIS矩阵的初始猜测如下
TRANS_GUESS=[.85.15;.1.9];
EMIS_GUESS=[.17.16.17.16.17.17;.6.08.08.08.0808];
下面做矩阵估计
[TRANS_EST2,EMIS_EST2]=hmmtrain(q,TRANS_GUESS,EMIS_GUESS)
TRANS_EST2=
0.22860.7714
0.00320.9968
EMIS_EST2=
0.14360.23480.18370.19630.23500.0066
0.43550.10890.11440.10820.11090.1220
Hmmtrain函数使用迭代算法调整TRANS_GUESS和EMIS_GUESS,使得基于调整的矩阵产
生更加类似于观察序列q的输出。当两次迭代中矩阵的变化不大时迭代结束。
如果如果算法在默认的最大迭代次数100次以内仍没结束,将返回迭代的最后结果并产
生一个警告。这时可以通过增加最大迭代次数来增大迭代的次数
hmmtrain(q,TRANS_GUESS,EMIS_GUESS,'maxiterations',maxiter)
maxiter是新设的最大迭代次数。
还可以改变默认的最小容忍值
hmmtrain(q,TRANS_GUESS,EMIS_GUESS,'tolerance',tol)
tol为新设的最小容忍值,增大它可以让算法更快的停止,但结果将不精确。
两个因素将降低hmmtrain函数估计结果的可靠性:
算法收敛到一个局部极大值,但此极大值并不是真正的状态转移矩阵与输出矩阵的位置。
对此,可以使用多个不同的初始猜测值来估计,使用最好的结果;
Seq序列太短,以至于不足以训练矩阵,对此,可以使用更长的q序列。
估计后验状态概率
输出序列q的后验状态概率是模型给出q序列时所处状态的后验概率。可以使用
hmmdecode函数计算后验状态概率:
PSTATES=hmmdecode(q,TRANS,EMIS)
随着序列长度的增加,序列的概率会趋向于0,所以一个足够长的序列的概率将小于计
算机的最小精度,为了避免此问题,hmmdecode还返回了概率的对数值。
改变初始状态分布
默认地,统计工具箱的隐马尔科夫模型函数是从状态1开始的。换句话说,初始状态的
分布集中在状态1处,其他位置为0。如果想给初始状态赋予一个概率分布p=[p1,p2,…,pM],
可以如下实现
1.创建一个M+1XM+1的辅助状态转移矩阵如下
其中的T是真正的状态转移矩阵,第一列包含M+1个0,向量p的和为1。
2.创建M+1XN的辅助输出矩阵如下
如果状态转移矩阵与输出矩阵式TRANS和EMIS,可以使用如下命令创建辅助矩阵
TRANS_HAT=[0p;zeros(size(TRANS,1),1)TRANS];
EMIS_HAT=[zeros(1,size(EMIS,2));EMIS];
本文发布于:2022-11-25 18:58:55,感谢您对本站的认可!
本文链接:http://www.wtabcd.cn/fanwen/fan/90/20278.html
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。
留言与评论(共有 0 条评论) |