HMM Package parameter estimation

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

HMM Package parameter estimation

Richard Philip
Hi,

I am having difficulties estimating the parameters of a HMM using the HMM
package. I have simulated a sequence of observations from a known HMM. When
I estimate the parameters of a HMM using these simulated observations the
parameters are not at all close to the known ones. I realise the estimated
parameters are not going to be exactly the same as the known/true
parameters, but these are nowhere close. Below is my code used. Any ideas
or possible suggestions regarding this issue would be greatly appreciated?


library(HMM)

## DECLARE PARAMETERS OF THE KNOWN MODEL
states = c(1,2,3)
symbols = c(1,2)
startProb = c(0.5,0.25,0.25)
transProb = matrix(c(0.8,0.05,0.15,0.2,0.6,0.2,0.2,0.3,0.5),3,3,TRUE)
emissionProb =  matrix(c(0.9,0.1,0.2,0.8,0.7,0.3), 3,2,TRUE)

# CREATE THE KNOWN MODEL
hmmTrue = initHMM(states, symbols, startProb, transProb , emissionProb)

# SIMULATE 1000 OBSERVATIONS OF THE KNOWN MODEL
observation = simHMM(hmmTrue, 1000)
obs = observation$observation

#ESTIMATE A MODEL USING THE OBSERVATIONS GENERATED FROM THE KNOWN MODEL
hmmInit = initHMM(states, symbols, c(1/3,1/3,1/3))
hmmFit = baumWelch(hmmInit, obs)


#The parameters of hmmTrue and hmmFit are not at all alike, why is this?



Kind Regards,
Richard

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: HMM Package parameter estimation

Rolf Turner-3

I think it's your starting values for the initial state probability
distribution,
i.e. c(1,1,1)/3 that cause the problem.  They seem to drop you into some
sort of local maximum/stationary point, a long way from the global maximum.

Try, e.g. c(4,2,1)/7; this gives me:

  hmmFit$hmm$emissionProbs
       symbols
states         1          2
      1 0.9385018 0.06149819
      2 0.7883591 0.21164092
      3 0.2279287 0.77207131

  hmmFit$hmm$transProbs
     to
from         1         2         3
    1 0.6925055 0.1239590 0.1835355
    2 0.2537700 0.5780679 0.1681621
    3 0.2455462 0.1190872 0.6353666

which look to be in "reasonable" agreement with the "true" values.
Note though that states 2 and 3 have been swapped.  This happens.

     cheers,

         Rolf Turner


On 16/04/13 13:13, Richard Philip wrote:

> Hi,
>
> I am having difficulties estimating the parameters of a HMM using the HMM
> package. I have simulated a sequence of observations from a known HMM. When
> I estimate the parameters of a HMM using these simulated observations the
> parameters are not at all close to the known ones. I realise the estimated
> parameters are not going to be exactly the same as the known/true
> parameters, but these are nowhere close. Below is my code used. Any ideas
> or possible suggestions regarding this issue would be greatly appreciated?
>
>
> library(HMM)
>
> ## DECLARE PARAMETERS OF THE KNOWN MODEL
> states = c(1,2,3)
> symbols = c(1,2)
> startProb = c(0.5,0.25,0.25)
> transProb = matrix(c(0.8,0.05,0.15,0.2,0.6,0.2,0.2,0.3,0.5),3,3,TRUE)
> emissionProb =  matrix(c(0.9,0.1,0.2,0.8,0.7,0.3), 3,2,TRUE)
>
> # CREATE THE KNOWN MODEL
> hmmTrue = initHMM(states, symbols, startProb, transProb , emissionProb)
>
> # SIMULATE 1000 OBSERVATIONS OF THE KNOWN MODEL
> observation = simHMM(hmmTrue, 1000)
> obs = observation$observation
>
> #ESTIMATE A MODEL USING THE OBSERVATIONS GENERATED FROM THE KNOWN MODEL
> hmmInit = initHMM(states, symbols, c(1/3,1/3,1/3))
> hmmFit = baumWelch(hmmInit, obs)
>
>
> #The parameters of hmmTrue and hmmFit are not at all alike, why is this?

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: HMM Package parameter estimation

Ingmar Visser
It seems that indeed providing other starting values initiates iterations
to take place. However, more worrisome is that the does not seem to
converge, even when upping the number of iterations. Below I run a 2 state
model on the same as well for comparison (I have added set.seed statements
to make the output exactly reproducible).
For 2 and 3 state models I get (see code below):

> fm2
Convergence info: Log likelihood converged to within tol. (relative change)
'log Lik.' -585.7628 (df=5)
AIC:  1181.526
BIC:  1206.064
> fm3
Convergence info: 'maxit' iterations reached in EM without convergence.
'log Lik.' -585.5807 (df=11)
AIC:  1193.161
BIC:  1247.147

hth, Ingmar

set.seed(11)
# SIMULATE 1000 OBSERVATIONS OF THE KNOWN MODEL
observation = simHMM(hmmTrue, 1000)
obs = observation$observation

#ESTIMATE A MODEL USING THE OBSERVATIONS GENERATED FROM THE KNOWN MODEL
hmmInit = initHMM(states, symbols, c(4,2,1)/7 )
# # hmmFit = baumWelch(hmmInit, obs)

hmmFit = baumWelch(hmmInit, obs, maxI=200)

library(depmixS4)

m2 <- depmix(obs~1,family=multinomial("identity"),ns=2,nt=1000)
set.seed(12)
fm2 <- fit(m2)

m3 <- depmix(obs~1,family=multinomial("identity"),ns=3,nt=1000)
set.seed(13)
fm3 <- fit(m3)



On Tue, Apr 16, 2013 at 11:53 AM, Rolf Turner <[hidden email]>wrote:

>
> I think it's your starting values for the initial state probability
> distribution,
> i.e. c(1,1,1)/3 that cause the problem.  They seem to drop you into some
> sort of local maximum/stationary point, a long way from the global maximum.
>
> Try, e.g. c(4,2,1)/7; this gives me:
>
>  hmmFit$hmm$emissionProbs
>       symbols
> states         1          2
>      1 0.9385018 0.06149819
>      2 0.7883591 0.21164092
>      3 0.2279287 0.77207131
>
>  hmmFit$hmm$transProbs
>     to
> from         1         2         3
>    1 0.6925055 0.1239590 0.1835355
>    2 0.2537700 0.5780679 0.1681621
>    3 0.2455462 0.1190872 0.6353666
>
> which look to be in "reasonable" agreement with the "true" values.
> Note though that states 2 and 3 have been swapped.  This happens.
>
>     cheers,
>
>         Rolf Turner
>
>
>
> On 16/04/13 13:13, Richard Philip wrote:
>
>> Hi,
>>
>> I am having difficulties estimating the parameters of a HMM using the HMM
>> package. I have simulated a sequence of observations from a known HMM.
>> When
>> I estimate the parameters of a HMM using these simulated observations the
>> parameters are not at all close to the known ones. I realise the estimated
>> parameters are not going to be exactly the same as the known/true
>> parameters, but these are nowhere close. Below is my code used. Any ideas
>> or possible suggestions regarding this issue would be greatly appreciated?
>>
>>
>> library(HMM)
>>
>> ## DECLARE PARAMETERS OF THE KNOWN MODEL
>> states = c(1,2,3)
>> symbols = c(1,2)
>> startProb = c(0.5,0.25,0.25)
>> transProb = matrix(c(0.8,0.05,0.15,0.2,0.**6,0.2,0.2,0.3,0.5),3,3,TRUE)
>> emissionProb =  matrix(c(0.9,0.1,0.2,0.8,0.7,**0.3), 3,2,TRUE)
>>
>> # CREATE THE KNOWN MODEL
>> hmmTrue = initHMM(states, symbols, startProb, transProb , emissionProb)
>>
>> # SIMULATE 1000 OBSERVATIONS OF THE KNOWN MODEL
>> observation = simHMM(hmmTrue, 1000)
>> obs = observation$observation
>>
>> #ESTIMATE A MODEL USING THE OBSERVATIONS GENERATED FROM THE KNOWN MODEL
>> hmmInit = initHMM(states, symbols, c(1/3,1/3,1/3))
>> hmmFit = baumWelch(hmmInit, obs)
>>
>>
>> #The parameters of hmmTrue and hmmFit are not at all alike, why is this?
>>
>
> ______________________________**________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/**listinfo/r-help<https://stat.ethz.ch/mailman/listinfo/r-help>
> PLEASE do read the posting guide http://www.R-project.org/**
> posting-guide.html <http://www.R-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
>

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.