by David Huen

We will now go through the source of one of the examples in the demos/dp directory: Dice.java (contributed by Samiul Hasan).

Introduction

The program implements the “occasionally dishonest casino” example used in the book “Biological Sequence Analysis” by R. Durbin, S. Eddy, A. Krogh, G. Mitchison.

Basically, it conceives a casino with two dice, one fair and one loaded. The fair die lands on any of its sides equal probability while the loaded die yields “6” half the time, all the other sides being of equal probability. These probabilities represent the emission distribution of the fair die state and the loaded die states respectively.

The casino switches between using the fair die and the loaded die periodically. When on the fair die, the probability that the next throw is with the fair die too is 0.95. Similarly, when on the loaded die, the probability of continuing with it is 0.90. These probabilities yield the transition distributions of the states.

The HMM as modelled in the code is slightly modified from the above description with the inclusion of a MagicalState. This state is used to represent the start and end of the states of the model. The transition from the MagicalState to the fair die state occurs with a probability of 0.8 while the transition to the loaded die state occurs with a probability of 0.2. A termination condition was also introduced to allow transitions from the fair die and loaded die states to the Magical state with a probability of 0.01.

Code

The core of the program is the createCasino() method. This creates an instance of the MarkovModel class that implements the model.

  public static MarkovModel createCasino() {
    Symbol[] rolls=new Symbol[6];

    //set up the dice alphabet
    SimpleAlphabet diceAlphabet=new SimpleAlphabet();
    diceAlphabet.setName("DiceAlphabet");

    for(int i=1;i<7;i++) {
      try {
        rolls[i-1]= AlphabetManager.createSymbol((char)('0'+i),""+i,Annotation.EMPTY_ANNOTATION);
        diceAlphabet.addSymbol(rolls[i-1]);
      } catch (Exception e) {
        throw new NestedError(
          e, "Can't create symbols to represent dice rolls"
        );
      }
    }

A Symbol array rolls is created to hold the symbols generated by AlphabetManager to represent the outcomes of the dice. An alphabet is also defined over these symbols.

Next, distributions representing the emission probabilities of the fair die and loaded die states are created (named fairD and loadedD respectively). The die states themselves are then created as SimpleEmissionStates, fairS and loadedS respectively.

You will observe an int array advance with a single value of 1. In a single-head HMM like ours, there is only one generated sequence and in our case, we progress along this sole sequence a single position per transition in the model. In multihead HMMs, there will be multiple sequences generated by the HMM and it is possible that the increment through the different sequences might be different. For example, single-stepping a protein sequence amounts to an increment of three on its corresponding DNA sequence.

int [] advance = { 1 };
Distribution fairD;
Distribution loadedD;
try {
  fairD = DistributionFactory.DEFAULT.createDistribution(diceAlphabet);
  loadedD = DistributionFactory.DEFAULT.createDistribution(diceAlphabet);
} catch (Exception e) {
  throw new NestedError(e, "Can't create distributions");
}
EmissionState fairS = new SimpleEmissionState("fair", Annotation.EMPTY_ANNOTATION, advance, fairD);
EmissionState loadedS = new SimpleEmissionState("loaded", Annotation.EMPTY_ANNOTATION, advance, loadedD);

The HMM is then created with these states:-

    SimpleMarkovModel casino = new SimpleMarkovModel(1, diceAlphabet, "Casino");
    try {
      casino.addState(fairS);
      casino.addState(loadedS);
    } catch (Exception e) {
      throw new NestedError(e, "Can't add states to model");
    }

Next, we need to model the transitions between the states. We do this like so:-

    try {
      casino.createTransition(casino.magicalState(),fairS);
      casino.createTransition(casino.magicalState(),loadedS);
      casino.createTransition(fairS,casino.magicalState());
      casino.createTransition(loadedS,casino.magicalState());
      casino.createTransition(fairS,loadedS);
      casino.createTransition(loadedS,fairS);
      casino.createTransition(fairS,fairS);
      casino.createTransition(loadedS,loadedS);
    } catch (Exception e) {
      throw new NestedError(e, "Can't create transitions");
    }

Note the presence of a MagicalState that is returned by casino.magicalState(). This is inherent to the SimpleMarkovModel class and does not need to be created by the user.

The emission distributions fairD and loadedD we set up earlier need to be initialised. We do that here.

    try {
      for(int i=0;i<rolls.length;i++)   {
        fairD.setWeight(rolls[i],1.0/6.0);
        loadedD.setWeight(rolls[i], 0.1);
      }
      loadedD.setWeight(rolls[5],0.5);
    } catch (Exception e) {
      throw new NestedError(e, "Can't set emission probabilities");
    }

We also need to initialise the transition distributions. Note how this is done: the transition distribution of each state is requested from the model with a getWeights() and then updated with the required values by calling the getWeight() method of that distribution. It is not necessary thereafter to call setWeights() to pass the distribution for a state back to the model. This may seem strange but it is done this way because model object may use unique Distribution classes that cannot be replaced by a generic Distribution class for greater internal efficiency. Every state in the model needs to have its own transition distribution initialised appropriately.

    //set up transition scores.
    try {
      Distribution dist;

      dist = casino.getWeights(casino.magicalState());
      dist.setWeight(fairS, 0.8);
      dist.setWeight(loadedS, 0.2);

      dist = casino.getWeights(fairS);
      dist.setWeight(loadedS,               0.04);
      dist.setWeight(fairS,                 0.95);
      dist.setWeight(casino.magicalState(), 0.01);

      dist = casino.getWeights(loadedS);
      dist.setWeight(fairS,                 0.09);
      dist.setWeight(loadedS,               0.90);
      dist.setWeight(casino.magicalState(), 0.01);
    } catch (Exception e) {
      throw new NestedError(e, "Can't set transition probabilities");
    }

Having completed constructing the MarkovModel, all that remains is to return it to the caller.

    return casino;

Using the MarkovModel

Having created the MarkovModel, we create the corresponding dynamic programming object:-

      DP dp=DPFactory.DEFAULT.createDP(casino);

Now, at last, we have something we can use! To generate a sequence of dice throws with this model, we do:-

      StatePath obs_rolls = dp.generate(300);

      SymbolList roll_sequence = obs_rolls.symbolListForLabel(StatePath.SEQUENCE);

The generate() method generates a path through the model that emits 300 symbols and we turn that path into a SymbolList with the second line.

At this point, it will be worthwhile digressing briefly on the StatePath object. This object embodies an Alignment of the sequences emitted by the DP object. In a multihead object, multiple aligned sequences will be emitted. In our case, only a single sequence is emitted and that is accessed with the label StatePath.SEQUENCE. That sequence turns out to a run of dice throws from our occasionally dishonest casino.

Next, we want to test one of the DP algorithms in the DP object. We want to process the roll_sequence SymbolList we have just generated and use the Viterbi method to predict which die each of the throws might have arisen from.

To do this, we create an array of SymbolLists with only roll_sequence in it - ours is a single-head HMM - and apply the Viterbi algorithm using the model probabilities (ScoreType.PROBABILITY) for the computation (you could have also applied the null-model or log-odds probabilities here). This will yield the state path that has most support from the model and that state path is the model’s prediction of which die a particular result came from.

      SymbolList[] res_array = {roll_sequence};
      StatePath v = dp.viterbi(res_array, ScoreType.PROBABILITY);

All that remains is to print out the generated sequence and the actual state path from which it came (ie. which die) and the HMM estimate of the state path, v, (which is which die a particular throw came from as estimated by the HMM).

      //print out obs_sequence, output, state symbols.
      for(int i = 1; i <= obs_rolls.length()/60; i++) {
        for(int j=i*60; j

The first two print statements print the generated sequence and the actual state path by accessing the obs_rolls StatePath<code> with the <code>StatePath.SEQUENCE and StatePath.STATES respectively. The predicted state path comes from the third print block which accesses the v StatePath.

The output then looks like this:-

544552213525245666363632432522253566166546666666533666543261
fffffffffffflllllllllllfffffffffffflllllllllllllllllllffffff
ffffffffffffffffffffffffffffffffffllllllllllllllllllllffffff

363546253252546524422555242223224344432423341365415551632161
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

144212242323456563652263346116214136666156616666566421456123
fffffflllfffffffffffffffffffffffflfllllllllllllllllfffffffff
fffffffffffffffffffffffffffffffffffllllllllllllllllfffffffff

346313546514332164351242356166641344615135266642261112465663
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

The top line is the sequence emitted by our HMM when we made it generate 300 throws. The next is the state from which the throw came (f-fair l-loaded, these are the first letters of the labels “fair” and “loaded” we used when creating the SimpleEmissionState objects that represent the dice). The last is similar but this time from the StatePath v that is the result of the Viterbi algorithm. The performance is pretty on this occasion but it can vary widely!