Hi! I'm planning to use HMMBase.jl for biological

Git Product

Light

Multiple sequence with different length about hmmbase.jl HOT 17 OPEN

sosuts commented on June 6, 2024

Multiple sequence with different length

from hmmbase.jl.

Comments (17)

sosuts commented on June 6, 2024 1

I restarted to write a code for multiple observations.
I think I can send PR soon.

from hmmbase.jl.

sosuts commented on June 6, 2024 1

Thanks to your original code!
I'll try fixing things step by step.

from hmmbase.jl.

cossio commented on June 6, 2024 1

Hello! I need this feature for some data I have, where also I have multiple time-series of different lengths.

What's the current status? I saw the linked PR got closed, but I'm not sure why?

Thanks!

from hmmbase.jl.

maxmouchet commented on June 6, 2024

Hi,

Parameter estimation with multiple sequences (of potentially different lengths), is not currently supported in HMMBase. It's something that I would like to implement (in mle.jl), but I don't know when.

I'm not sure if variable length sequences are supported by MS_HMMBase.

I'll keep this issue open as a reminder to try to implement this :-)

from hmmbase.jl.

sosuts commented on June 6, 2024

Thank you for your reply :)

At least I can implement the algorithm with a specific distribution, but I don’t know how to make it work with an arbitrary distributions...

from hmmbase.jl.

maxmouchet commented on June 6, 2024

In HMMBase I make use of the Distributions.jl package for the pdf (likelihood) and fit_mle (parameter estimation) methods. This way I can handle any observations distributions that implement these methods.
I only need to compute the "responsibilities" (assignments probabilities) for each observations and each states, and re-estimate the transition matrix.
See https://github.com/maxmouchet/HMMBase.jl/blob/master/src/mle.jl#L58-L68 where I delegate to fit_mle.

from hmmbase.jl.

sosuts commented on June 6, 2024

Hi again,
I'm now trying to integrate my code to this package.

BTW, I'm not sure what you mean by responsibilities of each observations.

from hmmbase.jl.

maxmouchet commented on June 6, 2024

By "responsibilities" I mean P(Z_t = i | Y, θ), where Z_t is the hidden state at time t, Y the observations, and θ the model parameters.

from hmmbase.jl.

sosuts commented on June 6, 2024

Thank you.
I have another question.
I was thinking that we scale β with c calculated from α.

HMMBase.jl/src/messages.jl

Lines 56 to 67 in f992852

 @inbounds for t in T-1:-1:1 

 for j in OneTo(K) 

 for i in OneTo(K) 

 β[t,j] += β[t+1,i] * A[j,i] * L[t+1,i] 

 end 

 c[t+1] += β[t,j] 

 end 

 for j in OneTo(K) 

 β[t,j] /= c[t+1] 

 end 

 end

Do I need to calculate another c in β?

from hmmbase.jl.

maxmouchet commented on June 6, 2024

You're right, it's possible to use the same scaling vector c for α and β.
It's not done, for now, to keep the code simple :-)

from hmmbase.jl.

sosuts commented on June 6, 2024

Sorry for asking many questions!
Why do we have to subtract m from loglikelihood?

HMMBase.jl/src/mle.jl

Lines 39 to 43 in 2ebf644

 m = vec_maximum(view(LL, t + 1, :)) 

 c = 0.0 

 for i in OneTo(K), j in OneTo(K) 

 ξ[t, i, j] = α[t, i] * A[i, j] * exp(LL[t+1, j] - m) * β[t+1, j]

HMMBase.jl/src/messages.jl

Lines 103 to 123 in 2ebf644

 m = vec_maximum(view(LL, 1, :)) 

 for j in OneTo(K) 

 α[1, j] = a[j] * exp(LL[1, j] - m) 

 c[1] += α[1, j] 

 end 

 for j in OneTo(K) 

 α[1, j] /= c[1] 

 end 

 c[1] *= exp(m) + eps() 

 @inbounds for t = 2:T 

 m = vec_maximum(view(LL, t, :)) 

 for j in OneTo(K) 

 for i in OneTo(K) 

 α[t, j] += α[t-1, i] * A[i, j] 

 end 

 α[t, j] *= exp(LL[t, j] - m)

from hmmbase.jl.

maxmouchet commented on June 6, 2024

No worries!

This is the log-sum-exp trick : https://en.wikipedia.org/wiki/LogSumExp
It prevents exp(LL[t,j]) from overflowing.

from hmmbase.jl.

maxmouchet commented on June 6, 2024

Nice!

I recently cleaned-up the code by removing the logl keyword and the methods that do not use the log-likelihoods.
Basically everything is done using the log-likelihoods now.

Feel free to open a PR, and I'll help you if there are merge conflicts.

from hmmbase.jl.

sosuts commented on June 6, 2024

Hi.
I changed my codes to adapt to your new api.
I implemented some new functions assuming 2 situations;

multiple observations with same length
multiple observations with different(random) length

I haven't finished writing codes for multivariate model in situation 2.
This notebook is an example.

Is it ok to open a pr?
To be honest, this is my first time using github so I'm not sure when to open it...

from hmmbase.jl.

maxmouchet commented on June 6, 2024

I think the correct URL is https://nbviewer.jupyter.org/github/SosUts/HMMBase.jl/blob/multiple_sequences/notebooks/multiple%20sequences.ipynb :)

This looks very nice!
Thank you for your work :)

You can open a PR now, and I'll review the code.
It is still possible to push new commits to your branch after the PR is opened, so there is no problem to make further modifications.

from hmmbase.jl.

sosuts commented on June 6, 2024

I opened it. Thank you as always!

Edit:
I forgot to consider about the tests.
Should I change the tests, or should I close the PR and change the codes?

from hmmbase.jl.

maxmouchet commented on June 6, 2024

No worries, you can keep the PR open!
Every commit that you add to your branch will be added to the PR automatically.

I'm a bit busy this week, so I'll try to have a look at the PR this week-end, or the next one.
In any case your code looks clean :)

from hmmbase.jl.

Related Issues (20)

Recommend Projects

React

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

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

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

An Open Source Machine Learning Framework for Everyone
Django

The Web framework for perfectionists with deadlines.
Laravel

A PHP framework for web artisans
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.
Visualization

Some thing interesting about visualization, use data art
Game

Some thing interesting about game, make everyone happy.

Recommend Org

Facebook

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

Open source projects and samples from Microsoft.
Google

Google ❤️ Open Source for everyone.
Alibaba

Alibaba Open Source for everyone
D3

Data-Driven Documents codes.
Tencent

China tencent open source team.

	@inbounds for t in T-1:-1:1
	for j in OneTo(K)
	for i in OneTo(K)
	β[t,j] += β[t+1,i] * A[j,i] * L[t+1,i]
	end
	c[t+1] += β[t,j]
	end

	for j in OneTo(K)
	β[t,j] /= c[t+1]
	end
	end

	m = vec_maximum(view(LL, t + 1, :))
	c = 0.0

	for i in OneTo(K), j in OneTo(K)
	ξ[t, i, j] = α[t, i] * A[i, j] * exp(LL[t+1, j] - m) * β[t+1, j]

	m = vec_maximum(view(LL, 1, :))

	for j in OneTo(K)
	α[1, j] = a[j] * exp(LL[1, j] - m)
	c[1] += α[1, j]
	end

	for j in OneTo(K)
	α[1, j] /= c[1]
	end

	c[1] *= exp(m) + eps()

	@inbounds for t = 2:T
	m = vec_maximum(view(LL, t, :))

	for j in OneTo(K)
	for i in OneTo(K)
	α[t, j] += α[t-1, i] * A[i, j]
	end
	α[t, j] *= exp(LL[t, j] - m)