toolbox是什么

更新时间:2022-11-25 18:58:55 阅读: 评论:0


2022年11月25日发(作者:广州自考)

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小时内删除。

上一篇:restrict
下一篇:fraud
相关文章
留言与评论(共有 0 条评论)
   
验证码:
Copyright ©2019-2022 Comsenz Inc.Powered by © 专利检索| 网站地图