分享好友 编程语言首页 频道列表

How to use the HMM toolbox (Matlab)

matlab  2023-02-09 10:230

一、离散输出的隐马尔科夫模型(DHMM,HMM with discrete outputs)

最大似然参数估计EM(Baum Welch算法)

The script dhmm_em_demo.m gives an example of how to learn an HMM with discrete outputs. Let there be Q=2 states and O=3 output symbols. We create random stochastic matrices as follows.

O = 3;

Q = 2;

prior0 = normalise(rand(Q,1));

transmat0 = mk_stochastic(rand(Q,Q));

obsmat0 = mk_stochastic(rand(Q,O)); 

Now we sample nex=20 sequences of length T=10 each from this model, to use as training data.

T=10;   %序列长度

nex=20;  %样本序列数目

data = dhmm_sample(prior0, transmat0, obsmat0, nex, T);      

Here data is 20x10. Now we make a random guess as to what the parameters are,

prior1 = normalise(rand(Q,1)); %初始状态概率

transmat1 = mk_stochastic(rand(Q,Q)); %初始状态转移矩阵

obsmat1 = mk_stochastic(rand(Q,O)); %初始观测状态到隐藏状态间的概率转移矩阵

and improve our guess using 5 iterations of EM...

[LL, prior2, transmat2, obsmat2] = dhmm_em(data, prior1, transmat1, obsmat1, 'max_iter', 5);

%prior2, transmat2, obsmat2 为训练好后的初始概率,状态转移矩阵及混合状态概率转移矩阵

LL(t) is the log-likelihood after iteration t, so we can plot the learning curve.

序列分类

To evaluate the log-likelihood of a trained model given test data, proceed as follows:

loglik = dhmm_logprob(data, prior2, transmat2, obsmat2)  %HMM测试

Note: the discrete alphabet is assumed to be {1, 2, ..., O}, where O = size(obsmat, 2). Hence data cannot contain any 0s.

To classify a sequence into one of k classes, train up k HMMs, one per class, and then compute the log-likelihood that each model gives to the test sequence; if the i'th model is the most likely, then declare the class of the sequence to be class i.

Computing the most probable sequence (Viterbi)

First you need to evaluate B(i,t) = P(y_t | Q_t=i) for all t,i:

B = multinomial_prob(data, obsmat);

Then you can use

[path] = viterbi_path(prior, transmat, B) 

 

二、具有高斯混合输出的隐马尔科夫模型(GHMM,HMM with mixture of Gaussians outputs)

Maximum likelihood parameter estimation using EM (Baum Welch)

Let us generate nex=50 vector-valued sequences of length T=50; each vector has size O=2.

O = 2;

T = 50;

nex = 50;

data = randn(O,T,nex);

Now let use fit a mixture of M=2 Gaussians for each of the Q=2 states using K-means.

M = 2;

Q = 2;

left_right = 0;

 

prior0 = normalise(rand(Q,1));

transmat0 = mk_stochastic(rand(Q,Q));

 

[mu0, Sigma0] = mixgauss_init(Q*M, reshape(data, [O T*nex]), cov_type);

mu0 = reshape(mu0, [O Q M]);

Sigma0 = reshape(Sigma0, [O O Q M]);

mixmat0 = mk_stochastic(rand(Q,M));

 

Finally, let us improve these parameter estimates using EM.

[LL, prior1, transmat1, mu1, Sigma1, mixmat1] = ...

    mhmm_em(data, prior0, transmat0, mu0, Sigma0, mixmat0, 'max_iter', 2);

Since EM only finds a local optimum, good initialisation is crucial. The initialisation procedure illustrated above is very crude, and is probably not adequate for real applications... Click here for a real-world example of EM with mixtures of Gaussians using BNT.

What to do if the log-likelihood becomes positive?

It is possible for p(x) > 1 if p(x) is a probability density function, such as a Gaussian. (The requirements for a density are p(x)>0 for all x and int_x p(x) = 1.) In practice this usually means your covariance is shrinking to a point/delta function, so you should increase the width of the prior (see below), or constrain the matrix to be spherical or diagonal, or clamp it to a large fixed constant (not learn it at all). It is also very helpful to ensure the components of the data vectors have small and comparable magnitudes (use e.g., KPMstats/standardize).

This is a well-known pathology of maximum likelihood estimation for Gaussian mixtures: the global optimum may place one mixture component on a single data point, and give it 0 covariance, and hence infinite likelihood. One usually relies on the fact that EM cannot find the global optimum to avoid such pathologies.

What to do if the log-likelihood decreases during EM?

Since I implicitly add a prior to every covariance matrix (see below), what increases is loglik + log(prior), but what I print is just loglik, which may occasionally decrease. This suggests that one of your mixture components is not getting enough data. Try a better initialization or fewer clusters (states).

What to do if the covariance matrix becomes singular?

Estimates of the covariance matrix often become singular if you have too little data, or if too few points are assigned to a cluster center due to a bad initialization of the means. In this case, you should constrain the covariance to be spherical or diagonal, or adjust the prior (see below), or try a better initialization.

How do I add a prior to the covariance matrix?

Buried inside of KPMstats/mixgauss_Mstep you will see that cov_prior is initialized to 0.01*I. This is added to the maximum likelihood estimate after every M step. To change this, you will need to modify the mhmm_em function so it calls mixgauss_Mstep with a different value.

Sequence classification

To classify a sequence (e.g., of speech) into one of k classes (e.g., the digits 0-9), proceed as in the DHMM case above, but use the following procedure to compute likelihood:

loglik = mhmm_logprob(data, prior, transmat, mu, Sigma, mixmat);

Computing the most probable sequence (Viterbi)

First you need to evaluate B(t,i) = P(y_t | Q_t=i) for all t,i:

B = mixgauss_prob(data(:,:,ex), mu, Sigma, mixmat);

where data(:,:,ex) is OxT where O is the size of the observation vector. Finally, use

[path] = viterbi_path(prior, transmat, B);

三、具有高斯输出的HMM

This is just like the mixture of Gaussians case, except we have M=1, and hence there is no mixing matrix.

Online EM for discrete HMMs/ POMDPs

For some applications (e.g., reinforcement learning/ adaptive control), it is necessary to learn a model online. The script dhmm_em_online_demo gives an example of how to do this.

 

查看更多关于【matlab】的文章

展开全文
相关推荐
反对 0
举报 0
评论 0
图文资讯
热门推荐
优选好物
更多热点专题
更多推荐文章
如何在Abaqus的python中调用Matlab程序
目录1. 确定版本信息2. 备份python3. 设置环境变量4. 安装程序5. 调试运行参考资料Abaqus2018操作系统Win10 64位Python版本2.7(路径C:\SIMULIA\CAE\2018\win_b64\tools\SMApy\python2.7)2. 备份python将上述的“python2.7”文件夹复制出来,避免因操作错误

0评论2023-03-16608

ROS与Matlab系列:一个简单的运动控制 基于matlab的运动控制系统
转自:http://blog.exbot.net/archives/2594Matlab拥有强大的数据处理、可视化绘图能力以及众多成熟的算法函数,非常适合算法开发;在控制系统设计中,Simulink也是普遍使用的设计和仿真工具。而ROS系统,则是一种新的标准化机器人系统软件框架。通过ROS,你

0评论2023-02-10920

matlab 遍历结构体struc的成员
MATLAB中专门用于对结构数组的操作的函数并不多,通过 help datatypes获取数据类型列表,可以看到其中的结构数据类型的有关的函数,主要如表4.3.1所示。表4.3.1 结构数组的操作函数函数名             功能描述 deal                 把输入处

0评论2023-02-09712

matlab编程如何换行 matlab怎么换行
空格+三个点+逗号

0评论2023-02-09691

C/C++中调用matlab引擎计算 matlab转c
原帖地址:http://blog.sina.com.cn/s/blog_6adcb3530101cvot.html一,在linux环境使用matlab引擎必须先进行一些必要的配置1,matlab引擎依赖/bin/csh启动,所以不管你使用何种shell,都必须安装csh。**2,matlab引擎依赖的动态库文件目录必须在系统当前的

0评论2023-02-09451

matlab几何纠正,间接法,双线性内插
超简洁,超级快,两个文件datapre.m文件代码:global X;global Y;global A;global l;global i;global I;global m;global n;global k;global dx;  A=[];l=[];i=0;m=[];n=[];dx=[]; fig=figure;subplot(1,2,1);I=imread('编程实习2-待纠正图像.bmp');imshow(I

0评论2023-02-09482

Matlab 根号的输入
二次根号:sqrt(a)或a^0.5三次根号:x^(1/3)或者x.^(1/3)根据x的数据结构类型矩阵、数组需要.^

0评论2023-02-09969

命令视频Matlab下查看摄像头设备信息
时间紧张,先记一笔,后续优化与完善。    应用如下2个命令:         info = imaqhwinfo('winvideo')    每日一道理 俄国作家契诃夫说:“有大狗,有小狗,小狗不该因为大狗的存在而心慌意乱。所有的狗都应该叫,就让他各自用上帝给他的声音

0评论2023-02-09957

permutation 随机置换检验的Matlab程序
假定a为某指标在10个样本中的值,5个一组,看以两组均值的差为例(统计量),随机置换检验程序 example: a: 230 -1350 -1580 -400 -760 970 110 -50 -190 -200v1=sum(a(1:5))/5;v2=sum(a(6:10))/5;T=abs(v1-v2);x=perms(a);      %矩阵a的全排列(随机全

0评论2023-02-09662

MATLAB学习1 之画图函数
ezplot适用条件“ezplot”命令可以用于显函数、隐函数和参数方程作图。不同函数的使用格式显函数y=f(x),ezplot函数的调用格式为ezplot(f, [xmin xmax]);              例:ezplot(\'sin(10*pi*x)/x\',[1 2]);%画出函数曲线隐函数f(x,y)=0,ezplot函数的

0评论2023-02-09591

Matlab 之 数据元素访问
Matlab的含义是矩阵实验室,其特征之一就是数据的向量化操作,借此提升软件运行效率。那么,必然会涉及数据元素的访问。Matlab主要支持下面一些形式的访问:(1)array-inde: A(i)(2)cell-index: C{i}(3)struct field: S.fieldname不同的访问方式,效

0评论2023-02-09994

matlab遍历文件夹下所有图片和遍历所有子文件夹下图片
做图像处理实验,经常需要遍历当前文件下所有图片。matlab当然很早就考虑了这个问题,库函数dir就是完成这个工作的。函数返回的是一个存放所有目录下文件信息的结构体,通过遍历结构体就可以达到访问所有文件的目的了。具体实现见下面程序:imgPath = 'E:/ima

0评论2023-02-09630

更多推荐