maxmouchet / hmmbase.jl Goto Github PK
View Code? Open in Web Editor NEWHidden Markov Models for Julia.
License: MIT License
Hidden Markov Models for Julia.
License: MIT License
Hi, thanks for the great package.
The current release for Distributions.jl is v0.25.37.
Is it possible to include this in the compat entry?
Thanks in advance.
fit_mle!
fit_mle
mle_step
messages*_log
, mle_step
and fit_mle!
Hi!
I'm planning to use HMMBase.jl
for biological usage.
Specifically, I want to estimate a parameters (e.g. speed) of an object which has multiple states.
My data have multiple coordinate data with different length.
I am new to julia so I'm not sure if HMMBase.jl
or MS_HMMBase.jl
supports this kind of analysis.
If not, are there any future plans?
Thank you in advance!
I'm sorry for a very silly question.
March 11: The results below are outdated/wrong, please see benchmark/ instead.
Performance @ 3856e94. Tested on 10k observations.
model = hmm.GaussianHMM(n_components=3, covariance_type="full")
model.startprob_ = np.array([.6, .3, .1])
model.transmat_ = np.array([[.7, .2, .1], [.3, .5, .2], [.3, .3, .4]])
model.means_ = np.array([[0.0, 0.0], [3.0, -3.0], [5.0, 10.0]])
model.covars_ = np.tile(np.identity(2), (3, 1, 1))
π0 = [.6, .3, .1]
P = [.7 .2 .1; .3 .5 .2; .3 .3 .4]
D = [MvNormal([0.0, 0.0], [1.0, 1.0]), MvNormal([3.0, -3.0], [1.0, 1.0]), MvNormal([5.0, 10.0], [1.0, 1.0])]
hmm = HMM(π0, P, D)
Times in ms. Mean time for Python via %timeit
, median for Julia via @benchmark
.
/ | hmmlearn | pyhsmm | HMMBase |
---|---|---|---|
Sample 10k | 981 | 189 | 1.818 |
lls | 1.11 | 3.562 | |
forward (w/o lls computation) | 5.08 | 1.510 | |
backward (w/o lls computation) | 5.02 | 1.402 | |
Viterbi (w/o lls computation) | 0.319 | 38.683 | |
Viterbi | 1.99 | 78.2 | 45.182 |
Just finished benchmarking what I hope is the final commit of MS_HMMBase, where I attempted to add multiple observations to the mle_step function and subfunctions, without compromising performance too much. Here are the results from running https://github.com/mmattocks/MS_HMMBase.jl/blob/master/test/benchmarks.jl:
[ Info: Judging log_likelihoods
BenchmarkTools.TrialJudgement:
time: +14.47% => regression (5.00% tolerance)
memory: -19.52% => improvement (1.00% tolerance)
[ Info: Judging messages_forwards_log
BenchmarkTools.TrialJudgement:
time: +1.04% => invariant (5.00% tolerance)
memory: +0.01% => invariant (1.00% tolerance)
[ Info: Judging messages_backwards_log
BenchmarkTools.TrialJudgement:
time: +1.35% => invariant (5.00% tolerance)
memory: +6.46% => regression (1.00% tolerance)
[ Info: Judging mle_step
BenchmarkTools.TrialJudgement:
time: -61.12% => improvement (5.00% tolerance)
memory: -29.02% => improvement (1.00% tolerance)
If you find any useful lines for HMMBase, please feel free to integrate them, obviously!
Also, if you have time, could you explain how messages_backwards_log works? It seems the views and the temp variable really make this much faster than simpler constructions might suggest- I haven't had much luck with views improving speed in other functions. Thanks very much for your time Max!!
Hi,
First of all, thank you for this nice package. I would like to use this package in a research project but I can not recover the original parameters when I use fit_mle. I create synthetic data using the function rand() and then I try to fit it with different initial conditions. For a HMM with bernulli distribution. The Transition matrix is correctly fit it but the fitted emission probabilities aare close to 1. Probably I am doing something wrong but I can not find my error. Here my code:
P=[0.7 0.3
0.1 0.9]
T=[0.7 0.3
0.1 0.9]
NDataSets=100
Nconditions=20
Nstates=2
Nout=2
Ntrials=1000
LlOriginal=zeros(NDataSets)
LlBest=zeros(NDataSets)
TBest=zeros(NDataSets,Nstates,Nstates)
PiBest=zeros(NDataSets,Nstates)
PBest=zeros(NDataSets,Nstates,Nout)
TFinal=zeros(Nconditions,Nstates,Nstates)
PiFinal=zeros(Nconditions,Nstates)
PFinal=zeros(Nconditions,Nstates,Nout)
Ll=zeros(Nconditions)
for idata in 1:NDataSets
println("iDataSet: ",idata)
#create a new hmm object
hmm = HMM(T[:,:], [Categorical(P[1,:]), Categorical(P[2,:])])
println(hmm.A)
println(hmm.B[1].p)
println(hmm.B[2].p)
choice,state=rand(hmm, Ntrials, seq = true) #create synthetic data
LlOriginal[idata]=-loglikelihood(hmm,choice) #LL of the original parameters
for icondition in 1:Nconditions
#randomize initial condition
aux=rand(4)
hmm.A[1,1]=aux[1]
hmm.A[1,2]=1-aux[1]
hmm.A[2,1]=aux[2]
hmm.A[2,2]=1-aux[2]
hmm.B[1].p[1]=aux[3]
hmm.B[1].p[2]=1-aux[3]
hmm.B[2].p[1]=aux[4]
hmm.B[2].p[2]=1-aux[4]
hmm2,history=fit_mle(hmm, choice) #fit
#println(hmm2.A)
TFinal[icondition,:,:]=hmm2.A
PFinal[icondition,1,:]=hmm2.B[1].p
PFinal[icondition,2,:]=hmm2.B[2].p
PiFinal[icondition,:]=hmm2.a
Ll[icondition]=-loglikelihood(hmm2,choice)
end
#best parameters fit of the Nconditions initial conditions
imin=findall(x->x==minimum(Ll),Ll)[1]
TBest[idata,:,:]=TFinal[imin,:,:]
PBest[idata,:,:]=PFinal[imin,:,:]
PiBest[idata,:,:]=PiFinal[imin,:,:]
LlBest[idata]=Ll[imin]
end
plot([T[1,1],T[1,1]],[1,3],"r-")
plot([T[2,2],T[2,2]],[1,3],"r-")
plot([P[1,1],P[1,1]],[3,5],"b--")
plot([P[2,2],P[2,2]],[3,5],"b--")
Hi,
do you have an example for a multivariate features Gaussian Mixture Model HMM, similar to GMMHMM from the hmmlearn Python package?
Extending the docs with such an example would help a lot.
Many thanks.
Hi Max!
I saw on the README that your package is looking for a new maintainer, and I feel like I have the ideal background to take care of it.
Right now I'm putting the finishing touches to my doctoral dissertation, but once September comes I should be more available. I'm excited to explore new features for a hypothetical 2.0 version, like compatibility with arbitrary emission distributions or transition structures. It would also be nice to work on compatibility with autodiff for gradient-based estimation.
What do you say?
Hi,
thanks for this nice package!
I was wandering what is the advantage of having the subtype VariateForm
AbstractHMM{F<:VariateForm}
?
I feel there is no real necessity of having this since you can get it from the vector of distribution.
I would like to solve exercise c) and d) of the attached problem with HMMBase.jl. I'm quite new to HMM but as I understand, HMMBase needs some observation values to perform the forward, backward and viterbi-algorithm whereas the problem only provides probabilities. Therefore it seems I can't solve this problem with HMMBase. Please let me know if I misunderstand the problem and/or the usage of HMMBase.jl
I was wondering if it was intentional to have the kmeans initialization use fit_mle instead of an estimator a user may pass as an argument? Thanks!
Dict/json/...
I'm trying to fit an hmm using the MvLogNormal
distribution, but it returns the following:
A = randtransmat(2)
B = [MvLogNormal(μ_pos, Σ),
MvLogNormal(μ_neg, Σ)]
hmm = HMM(A, B)
size(hmm)
hmm, hist = fit_mle(hmm, annual_data)
ERROR: suffstats is not implemented for (MvLogNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}}, Array{Float64,2}, Array{Float64,1}).
If I change it to MvNormal
it works perfectly well.
Hi guys,
In your mle.jl file, on line 14 in the mle_step function, you have:
log_likelihoods = HMMBase.log_likelihoods(hmm, observations)
If you define the loglikelihoods function for both, univariate and multivariate, with a fixed type like:
function log_likelihoods(hmm::AbstractHMM{Univariate}, observations)::Array{Float64,2}
hcat(map(d -> logpdf.(d, observations), hmm.D)...)
end
you should archieve a drastic speed up. I had a similar problem and it did the trick. I am no expert, but I think the 'log_likelihoods' functions within the 'mle' function is not type stable.
Also, I have to say that your module is really well structured and intuitively written!
I'm coming in relatively fresh to HMMs and I'm having trouble matching up terms I'm seeing in other work.
"Emission matrix" doesn't appear in the documentation - I gather it can be provided in place of B
, but the documentation doesn't show that or describe its form.
The form the "Transition matrix" should take isn't explicitly documented.
I'd be happy to make a contribution here, once I understand them.
I'm trying to fit a HMM to my data and I've reduce the problem to a MWE, I'll update if I can reduce it further.
using HMMBase #v1.0.6
using Distributions
using Random
model = MixtureModel([Normal(0.0, 0.0), LogNormal(7.9, 0.49)], [0.65, 0.35])
Random.seed!(123)
data=rand(model, 10000)
# adapted from https://github.com/maxmouchet/HMMBase.jl/issues/25
# to fix suffstats error
function Distributions.fit_mle(D::Type{LogNormal{T}}, x::AbstractMatrix, w::AbstractVector) where {T}
LogNormal(fit_mle(Normal{T}, log.(x), w))
end
hmm = HMM([0.65 0.35; 0.35 0.65], [
Normal(0.0, 1.0),
LogNormal(7.9, 0.49)
])
fit_mle(hmm, rand(model, 10, 241), display=:iter)
Gives me the following error
Iteration 0: logtot = -46.973035945669366
ArgumentError: Normal: the condition σ >= zero(σ) is not satisfied.
Stacktrace:
[1] macro expansion at /home/tlnagy/.julia/packages/Distributions/jFoHB/src/utils.jl:6 [inlined]
[2] #Normal#101 at /home/tlnagy/.julia/packages/Distributions/jFoHB/src/univariate/continuous/normal.jl:37 [inlined]
[3] Normal at /home/tlnagy/.julia/packages/Distributions/jFoHB/src/univariate/continuous/normal.jl:37 [inlined]
[4] fit_mle at /home/tlnagy/.julia/packages/Distributions/jFoHB/src/univariate/continuous/normal.jl:353 [inlined]
[5] fit_mle(::Type{Normal{Float64}}, ::Array{Float64,2}, ::Array{Float64,1}; mu::Float64, sigma::Float64) at /home/tlnagy/.julia/packages/Distributions/jFoHB/src/univariate/continuous/normal.jl:380
[6] fit_mle(::Type{Normal{Float64}}, ::Array{Float64,2}, ::Array{Float64,1}) at /home/tlnagy/.julia/packages/Distributions/jFoHB/src/univariate/continuous/normal.jl:378
[7] update_B!(::Array{Distribution{Univariate,S} where S<:ValueSupport,1}, ::Array{Float64,2}, ::Array{Float64,2}, ::typeof(fit_mle)) at /home/tlnagy/.julia/packages/HMMBase/VeHho/src/mle.jl:77
[8] fit_mle!(::HMM{Univariate,Float64}, ::Array{Float64,2}; display::Symbol, maxiter::Int64, tol::Float64, robust::Bool, estimator::Function) at /home/tlnagy/.julia/packages/HMMBase/VeHho/src/mle.jl:118
[9] #fit_mle#6 at /home/tlnagy/.julia/packages/HMMBase/VeHho/src/mle_api.jl:22 [inlined]
[10] top-level scope at In[2]:23
[11] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091
Any thoughts? Thanks!
Hi,
Thank you for the nice work of creating the package.
I am trying to use the package to first learn the parameters by fit_mle. Then I used the Viterbi function with the learned model and the observations to predict the hidden states. So far I've been able to get the hmm model parameters by following code
k = 5
ob_dim = size(zmat, 2)
A = randtransmat(k)
∑ = [MvNormal(rand(ob_dim), ones(ob_dim)) for i = 1:k]
hmm = HMM(A, ∑)
hmm, history = fit_mle(hmm, zmat, display=:final, maxiter=200, init=:kmeans)
I first try to decode the observations
viterbi(hmm, zmat)
I am having the error "ERROR: BoundsError: attempt to access 11823×5 Array{Int64,2} at index [11822, 0]", which 11823 is the number of observations (size(zmat,1) or T).
Then I try to sample another pseudo observations by rand and then Viterbi like in the package document
z, y = rand(hmm, 500, seq = true) #ok
viterbi(hmm, y)
"ERROR: BoundsError: attempt to access 500×5 Array{Int64,2} at index [499, 0]". I got the "same" err again. I am not sure what's going wrong in my code. I really appreciate your help!
Thanks
Nevermind, I misread the definition of a variable! Apologies for the spam.
Hi,
I am new to EM algorithm for hmm. I want to use fit_mle
for a product distribution of Bernoulli variable (and Categorical). I can define the hmm, I can sample from it. When I try fit_mle
ran into that message:
ERROR: suffstats is not implemented for (Product{Discrete,Bernoulli{Float64},Array{Bernoulli{Float64},1}}, Array{Int64,2}, Array{Float64,1})
.
Do you think it would be hard to actually add that? Do you have indications?
Array{Float64,2}
)Implementations of forward/backward/MLE dependent on the transition matrix type.
dellison:
The reason that I’m asking is that in natural language processing, HMMs have been used for sequence labeling tasks like part-of-speech tagging. In this task, for example, a second-order HMM that conditions on the previous two tags generally seems to do better than a first-order model that just conditions on one previous tag.
maxmouchet:
I’m not very familiar with higher-order models, but right now I see two ways of implementing them:
Writing specialized versions of the functions (e.g. messages_forwards_2); Implementing more generic algorithms (such as belief propagation) to handle models of arbitrary orders in a single function;
The second option seems cleaner but I’m worried about the performance implications. On the other hand, macros could be used to generate specialized function for an arbitrary order instead.
As for the types, I can add a new parameter like AbstractHMM{F<:VariateForm, N} where N represents the order (1 by default).
Array(HMM), Dict, JSON, ...
c
and to compute unscaled messages (set c
to 1?)Hello,
Thank you for this very nice package. I am just learning HMM and I would like to know if is it possible to estimate a HMM with 2 hidden states for 3 observations like this:
hmm = HMM([0.9 0.1; 0.1 0.9], [MvNormal(ones(3)), MvNormal(ones(3))])
nstates, obs_dim = size(hmm)
y = rand(hmm, 1000)
I tried it but it gives an error of dimension mismatch:
ERROR: DimensionMismatch("tried to assign 3-element array to 1×2 destination")
Am I doing something wrong here?
Hi Max,
I am doing something a bit strange with a fork of this library (I need to learn many HMMs on very large multisequence datasets) here: https://github.com/mmattocks/MS_HMMBase.jl
I am presently writing tests to ensure my modified functions are equivalent to yours. In doing so, I noticed that the following lines in mle.jl appear to have a small indexing error:
new_π = sum(ξ, dims=1)[1,:,:]
new_π ./= sum(new_π, dims=2)
Rabiner 89 has the summation over t to T-1 for xi / gamma, so I think we want ξ[1:end-1,:,:] for the first line. Please let me know if this is incorrect!
Thanks for your work on this library by the way, I've learned a lot from it!
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.