Hidden Markov Model segmentation of different proportions of binary data

Question

I need to segment a sequence of 0s and 1s by their proportion at relatively large scales. As an example, let's define 5 different states that represent 5 different ratios of 1s & 0s.

Alphabet: 1 and 0

State      Definition                  emission prob.
state 0: 100% zeroes and 0% ones      0:0.999   1: 0.001
state 1: 75% zeroes and 25% ones      0:0.75     1: 0.25
state 2: 50% zeroes and 50% ones      0: 0.5     1: 0.5
state 3: 25% zeroes and 75% ones      0: 0.25    1: 0.25
state 4: 0% zeroes and 100% ones      0: 0.001   1: 0.999

Attempts: With all the transition probabilities that I've tried so far and the emissions of each state, the output of my model sequences of states is only either state 0 or state 4.

Example:

data (binary):

00000000000001111111111110000000000101010101010101010000001000100010011001000010

output I get no matter how I change the transition probs. (in states):

00000000000004444444444440000000000404040404040404040000004000400040044004000040

output I need (in states):

00000000000004444444444440000000000222222222222222220000001111111111111111111111

I have the impression that I am missing some basic theory rather than an implementation problem. For instance, I smoothed the data by aggregating obtaining the ratio of 1 vs 0 in a arbitrarily defined window, and in this way I can see the intermediate states between state 0 and state 4. Nevertheless, I don't want to smooth the real data as I need to justify then the smoothing window size.

Is using HMMs a good solution for this problem?

Answer 1

My response is in two parts. First, by changing the input (initial) transition probabilities, you can get something similar to what you'd like. Here's some R code demonstrating this for your example:

library(HMM)

States <- c("0","1","2","3","4")
Symbols <- c("0","1")
startProbs <- rep(0.2,5)
emissionProbs <- matrix(c(0.999,0.75,0.5,0.25,0.001,0.001,0.25,0.5,0.75,0.999),5,2)
transProbs <- matrix(0.025,5,5)
diag(transProbs) <- 0.9

hmm <- initHMM(States, Symbols, startProbs, transProbs, emissionProbs)
> print(hmm)
$States
[1] "0" "1" "2" "3" "4"

$Symbols
[1] "0" "1"

$startProbs
  0   1   2   3   4 
0.2 0.2 0.2 0.2 0.2 

$transProbs
    to
from     0     1     2     3     4
   0 0.900 0.025 0.025 0.025 0.025
   1 0.025 0.900 0.025 0.025 0.025
   2 0.025 0.025 0.900 0.025 0.025
   3 0.025 0.025 0.025 0.900 0.025
   4 0.025 0.025 0.025 0.025 0.900

$emissionProbs
      symbols
states     0     1
     0 0.999 0.001
     1 0.750 0.250
     2 0.500 0.500
     3 0.250 0.750
     4 0.001 0.999

With this initial transition matrix, we get the following probabilities for observations 8, 20, 30, and 40, which are in the middle (roughly) of sequences of 0, 1, 0, and 0,1,0,1... respectively:

obs <- as.character(c(0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,
         0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,0,0,1,0,
         0,0,1,0,0,0,1,0,0,1,1,0,0,1,0,0,0,0,1,0))

post <- posterior(hmm, obs)
> post[,c(8,20,30,40)]
      index
states           8          20          30         40
     0 0.934764162 0.000001395 0.725475508 0.00004174
     1 0.059970011 0.000724501 0.244379742 0.31189082
     2 0.004632750 0.006383026 0.028836815 0.56445433
     3 0.000631774 0.082112681 0.001305885 0.11840354
     4 0.000001303 0.910778397 0.000002049 0.00520957

As you can see, the max. probability states are 0, 4, 0, and 2 respectively, as you wished.

It may also help you out if you don't pick such extreme probabilities for states 0 and 4, perhaps choosing 0.95 / 0.05 instead of 0.999 / 0.001. This will make it easier to have higher transition probabilities out of a given state without winding up in states 0 and 4 all the time.

If you are considering alternatives to HMMs, you might consider a continuous state space model, which can be formulated as a generalized additive model. Using the mgcv package in R, this could be set up as follows:

library(mgcv)

obs <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,
                      0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,0,0,1,0,
                      0,0,1,0,0,0,1,0,0,1,1,0,0,1,0,0,0,0,1,0)
time <- seq(1,length(obs))

foo <- gam(obs~s(time),family="binomial")

> predict(foo,type="response")[c(8,20,30,40)]
                   8                   20                   30 
0.000000000000000222 0.999999999999999778 0.000277113887323986 
                  40 
0.540166858432701846 

As you can see, the probabilities line up pretty well with what you'd like. Obviously some tuning of the parameters in the smoothing term would likely be desirable.