CODEWITHSUNDEEP
simulationlinear-algebrahidden-markov-models
mac389
mac389

Reputation: 3133

Generating a set of emissions given a transition matrix and starting state in a hidden markov model

I have the transition matrix, emission matrix and starting state for a hidden Markov model. I want to generate a sequence of observations (emissions). However, I'm stuck on one thing.

I understand how to choose among two states (or emissions). If Event A probability x, then Event B (or, really not-A) occurs with probability 1-x. To generate a sequence of A's and B's, with a random number, rand, you do the following.

 for iteration in iterations:
      observation[iteration] <- A if rand < x else B

I don't understand how to extend this to more than two variables. For example, if three events occur such that Event A occurs with probability x1, Event B with x2 and Event C with 1-(x1+x2), then how do I extend the above pseudocode?

I didn't find the answer Googling. In fact I get the impression that I'm missing a basic fact that many of the notes online assume. :-/

Upvotes: 0

Views: 1922

Answers (2)

Zhubarb
Zhubarb

Reputation: 11885

The two value case is a binomial distribution, and you generate random draws from a binomial distribution (a series of coin flips, essentially).

For more than 2 variables, you need to draw samples from a multinomial distribution, which is simply a generalisation of the binomial distribution for n>2.

Regardless of what language you use, there should most likely be built-in functions to accomplish this task. Below is some code in Python, which simulates a set of observations and states given your hmm model object:

import numpy as np
def random_MN_draw(n, probs):
    """ get X random draws from the multinomial distribution whose probability is given by 'probs' """
    mn_draw = np.random.multinomial(n,probs) # do 1 multinomial experiment with the given probs with probs= [0.5,0.5], this is a coin-flip
    return np.where(mn_draw == 1)[0][0] # get the index of the state drawn e.g. 0, 1, etc.

def simulate(self, nSteps):
    """ given an HMM = (A, B1, B2 pi), simulate state and observation sequences """
    lenB= len(self.emission)
    observations = np.zeros((lenB, nSteps), dtype=np.int) # array of zeros
    states = np.zeros(nSteps)
    states[0] = self.random_MN_draw(1, self.priors) # appoint the first state from the prior dist
    for i in range(0,lenB): # initialise observations[i,0] for all observerd variables
        observations[i,0] = self.random_MN_draw(1, self.emission[i][states[0],:]) #ith variable array, states[0]th row

    for t in range(1,nSteps): # loop through t
        states[t] = self.random_MN_draw(1, self.transition[states[t-1],:]) # given prev state (t-1) pick what row of the A matrix to use

        for i in range(0,lenB): # loop through the observed variable for each t
            observations[i,t] = self.random_MN_draw(1, self.emission[i][states[t],:]) # given current state t, pick what row of the B matrix to use

    return observations,states

In pretty much every language, you can find equivalents of

np.random.multinomial()

for multinomial and other discrete or continuous distributions as built-in functions.

Upvotes: 1

dmuir
dmuir

Reputation: 644

One way would be 

 x<-rand()
  if x < x1 observation is A
  else if x < x1 + x2 observation is B
  else observation is C

Of course if you have a large number of alternatives it might be better to build a cumulative probability table (holding x1, x1+x2, x1+x2+x3 ...) and then do a binary search in that table given the random number. If you are willing to do more preprocessing, there is an even more efficient way, see here for example

Upvotes: 1

Related Questions