Hidden Markov Model algorithms

HMM.cpp

//
//  main.cpp
//  HMM
//
//  Created by Peter Goldsborough on 01/07/14.
//  Copyright (c) 2014 Peter Goldsborough. All rights reserved.
//

#include <iostream>
#include <vector>

typedef std::vector<double> Array;
typedef std::vector<std::vector<double>> Matrix;
typedef std::vector<int> observSeq;

// Bayes Theorem
void normalize(Array& arr)
{
    double a,b;
    
    a = arr[0] / (arr[0] + arr[1]);
    b = arr[1] / (arr[0] + arr[1]);
    
    arr[0] = a;
    arr[1] = b;
}

Matrix forward(Array stateM, const Matrix& transM,const Matrix& observM, const observSeq& obs)
{
    Matrix ret;
    
    ret.push_back(stateM);
    
    for(unsigned long o = 0; o < obs.size(); ++o)
    {
        Array temp;
        
        // old state probabilities * transition probabilities
        for(unsigned long col = 0, colEnd = transM[0].size(); col < colEnd; ++col)
        {
            double val = 0;
            
            for (unsigned long row = 0; row < stateM.size(); ++row)
            {
                val += stateM[row] * transM[row][col];
            }
            
            temp.push_back(val);
        }

        stateM = temp;
        
        // new state probabilities * observation probabilities
        for (unsigned long state = 0; state < stateM.size(); ++state)
        {
            stateM[state] *= observM[state][obs[o]];
        }
        
        normalize(stateM);
        
        ret.push_back(stateM);
    }
    
    return ret;
}

Matrix backward(const Matrix& transM,const Matrix& observM, const observSeq& obs)
{
    // Assume uniform probability first
    Array stateM(observM.size(),1);
    
    Matrix ret;
    
    ret.push_back(stateM);
    
    for(long o = obs.size() - 1; o >= 0; --o)
    {
        // We need to modify the transition values but not the originals
        Matrix tempM(transM);
        
        // Same goes for the state matrix
        Array tempA;
        
        // First we multiply the transition probabilites by the observation
        // probabilites (given the current observation 'o' in the sequence)
        for (unsigned long state = 0; state < observM.size(); ++state)
        {
            for (unsigned long col = 0; col < transM[0].size(); ++col)
            {
                tempM[state][col] *= observM[col][obs[o]];
            }
        }
        
        // Then multiply those new probabilites times the by the previous
        // state's probabilites, giving the current latent variable's
        // probabilities
        for(unsigned long row = 0; row < tempM.size(); ++row)
        {
            double val = 0;
            
            for (unsigned long col = 0; col < stateM.size(); ++col)
            {
                val += stateM[col] * tempM[row][col];
            }
            
            tempA.push_back(val);
        }
        
        stateM = tempA;
        
        // normalize via bayes
        normalize(stateM);
        
        ret.push_back(stateM);
    }
    
    // get actual order (from 0 to n)
    std::reverse(ret.begin(), ret.end());
    
    return ret;
}

Matrix forward_backward(Array stateM, const Matrix& transM,const Matrix& observM, const observSeq& obs)
{
    Matrix alphas = forward(stateM, transM, observM, obs);
    
    Matrix betas = backward(transM, observM, obs);
    
    // multiply the matrices
    for (unsigned long i = 0; i < alphas.size(); ++i)
    {
        for (unsigned long j = 0; j < stateM.size(); ++j)
        {
            alphas[i][j] *= betas[i][j];
        }
        
        // normalize the values again
        normalize(alphas[i]);
    }
    
    return alphas;
}

observSeq viterbi (Array stateM, const Matrix& transM, const Matrix& observM, observSeq obs)
{
    std::vector<observSeq> path;
    path.resize(obs.size());
    
    Matrix probab;
    probab.resize(obs.size());
    
    for(Matrix::iterator itr = probab.begin(), end = probab.end();
        itr != end;
        ++itr)
    { itr->resize(stateM.size()); }
    
    // initialize values
    
    for (unsigned short s = 0; s < stateM.size(); ++s)
    {
        probab[0][s] = stateM[s] * observM[s][obs[0]];
        
        observSeq o = {s};
        path[s] = o;
    }
    
    // for every observation given
    for (unsigned long t = 1; t < obs.size(); ++t)
    {
        // store temporarily as to not change path while
        // using it
        std::vector<observSeq> temppath;
        temppath.resize(obs.size());
        
        // for every state
        for(unsigned short s = 0; s < stateM.size(); ++s)
        {
            double p = 0; // best probability
            unsigned short state = 0; // best state
            
            // calculate all possible paths via the forward algorithm
            // previous state probabilities * transition probabiliites
            // * the observation probabilities. Store only the best.
            for(unsigned short i = 0; i < stateM.size(); ++i)
            {
                double val = probab[t-1][i] * transM[i][s] * observM[s][obs[t]];
                
                // update maximum
                if (val > p)
                {
                    p = val;
                    state = i;
                }
            }
            
            probab[t][s] = p;
            
            temppath[s] = path[state];
            temppath[s].push_back(s);
        }
        
        path = temppath;
    }
    
    // If there is only one item, it is the maximum
    unsigned short n = (obs.size() > 1) ? (obs.size() - 1) : 0;
    
    double p = 0;
    unsigned short state = 0;
    
    // return the best
    for(unsigned short i = 0; i < stateM.size(); ++i)
    {
        if (probab[n][i] > p)
        {
            p = probab[n][i];
            state = i;
        }
    }
    
    std::cout << p << std::endl;
    
    return path[state];
}

int main(int argc, char * argv[])
{
    Matrix transM = {{0.7,0.3},{0.3,0.7}};
    Matrix observM = {{0.9,0.1},{0.2,0.8}};
    
    Array stateM = {0.5,0.5};
    
    observSeq obs = {0,0,1,0,0};
    
    Matrix m = forward_backward(stateM,transM, observM, obs);
    
    for (Matrix::const_iterator itr = m.begin(), end = m.end();
         itr != end;
         ++itr)
    {
        std::cout << itr - m.begin() << ": ";
        
        for (Array::const_iterator jtr = itr->begin(), jend = itr->end();
             jtr != jend;
             ++jtr)
        {
            std::cout << *jtr << "\t";
        }
        
        std::cout << std::endl;
    }
    
}