Git Product home page Git Product logo

hmmbase.jl's Issues

HMMBase v1.0

  • Implement the MLE estimator (Distributions::fit_mle! for the M step)
    • Implement fit_mle!
    • Implement fit_mle
    • Handle multivariate distributions in mle_step
    • Test messages*_log, mle_step and fit_mle!
  • Cleanup the Viterbi implementation (variable names, ...)
  • Benchmark against other implementations (C, Python, hmmlearn, ...) - see #2
  • Better tests (in particular check that Viterbi output is correct)
  • Compute stationnary distribution (if it exists)
  • EM: K-means initialization

Multiple sequence with different length

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.

Performance

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

Benchmarks vs MS_HMMBase & question about messages_backwards_log

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!!

I can not recover the original parameters

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--")

Application for new maintainer

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?

AbstractHMM subtyping

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.

HMM with observations as probabilities

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

hmm_example.pdf

Unable to fit using Multivariate LogNormal distribution

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.

Em Algorithm Speedup

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!

Improve documentation for novices

I'm coming in relatively fresh to HMMs and I'm having trouble matching up terms I'm seeing in other work.

  1. "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.

  2. The form the "Transition matrix" should take isn't explicitly documented.

I'd be happy to make a contribution here, once I understand them.

Running into an argument error in Distributions when using fit_mle

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!

viterbi(hmm, y) got "ERROR: BoundsError: attempt to access T×5 Array{Int64,2} at index [T-1, 0]"

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

Product of Discrete, Bernoulli

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?

Ideas

Specialized transition matrix types

  • Classical, dense, transition matrix (Array{Float64,2})
  • Left-right
  • Sparse

Implementations of forward/backward/MLE dependent on the transition matrix type.

Higher-order models

As per https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/2:

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).

Implement serialization/deserialization

Array(HMM), Dict, JSON, ...

Others

2 state and 3 observation HMM

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?

Error in mle.jl?

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!

TagBot trigger issue

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!

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.