classification/LeftToRightHMM.cpp

Go to the documentation of this file.

#include "LeftToRightHMM.h"

LeftToRightHMM::LeftToRightHMM(size_t nStates, size_t alphabetSize) :
    m_nStates(nStates),
    m_transitionProps(nStates, nStates),
    m_outputProps(nStates, alphabetSize),
    m_trelis(1,1)
{
    // init transitions and outputs with equally distributed values
    size_t i, j;
    for (i = 0; i < nStates; ++i) {
        for (j = i; j < nStates; ++j) {
            m_transitionProps[i][j] = 1.0 / static_cast<double>(nStates - i);
        }
        for (j = 0; j < alphabetSize; ++j) {
            m_outputProps[i][j] = 1 / static_cast<double>(alphabetSize);
        }
    }
}

double LeftToRightHMM::evaluate(Observation &obs)
{
    forwardTrelis(obs);
    return m_trelis[obs.size()][m_nStates-1];
}

void LeftToRightHMM::train(Observation &obs, size_t maxIterations)
{
    Matrix<double> alpha(obs.size() + 1, m_nStates, 0);
    Matrix<double> beta(obs.size() + 1, m_nStates, 0);
    
    std::vector<Matrix<double> > xi;
    for (size_t i = 0; i < obs.size(); ++i) {
        xi.push_back(Matrix<double>(m_nStates, m_nStates, 0));
    }

    for (size_t iteration = 0; iteration < maxIterations; ++iteration) {
        calcXi(obs, alpha, beta, xi);
        adjustParameters(obs, xi);
    }
}

void LeftToRightHMM::calcXi(Observation &obs, 
        Matrix<double> &alpha, 
        Matrix<double> &beta, 
        std::vector<Matrix<double> > &xi)
{
    Observation::iterator it;
    size_t t, i, j;
    double obsPropability = evaluate(obs);
    // forward trelis is still stored in m_trelis - no need to calc a second time
    alpha = m_trelis;
    backwardTrelis(obs);
    beta = m_trelis;
    for (t = 0, it = obs.begin(); t < obs.size() && it != obs.end(); ++t, ++it) {
        for (i = 0; i < m_nStates; ++i) {
            for (j = 0; j < m_nStates; ++j) {
                xi[t][i][j] = alpha[t][i] 
                            * m_transitionProps[i][j]
                            * m_outputProps[j][*it]
                            * beta[t+1][j]
                            / obsPropability;
            }
        }
    }
}

void LeftToRightHMM::adjustParameters(Observation &obs,
        std::vector<Matrix<double> > &xi)
{
    Observation::iterator it;
    size_t t, i, j, k;
    double divident, divisor;
    // adjust transition propabilites
    for (i = 0; i < m_nStates; ++i) {
        for (j = 0; j < m_nStates; ++j) {
            divident = divisor = 0;
            for (t = 0; t < obs.size(); ++t) {
                divident += xi[t][i][j];
                for (k = 0; k < m_nStates; ++k) {
                    divisor += xi[t][i][k];
                }
                m_transitionProps[i][j] = divident / divisor;
            }
        }
    }
    
    // adjust output propabilites
    for (j = 0; j < m_nStates; ++j) {
        for (k = 0; k < m_outputProps.getM(); ++k) {
            divident = divisor = 0;
            for (i = 0; i < m_nStates; ++i) {
                for (t = 0; t < obs.size(); ++t) {
                    if (obs[t] == static_cast<int>(k))
                        divident += xi[t][i][j];
                    divisor += xi[t][i][j];
                }
            }
            m_outputProps[j][k] = divident / divisor;
        }
    }
}

Matrix<double> &LeftToRightHMM::forwardTrelis(Observation &obs)
{
    m_trelis = Matrix<double>(obs.size() + 1, m_nStates, 0);
    // only state 1 (= start state) has propability 1
    m_trelis[0][0] = 1;

    // alpha_t(j) = (sum{i=0}{N}{alpha_{t-1}(i) * a_{i,j}}) * b_j(O_t)
    Observation::iterator it;
    size_t t = 1, j, i;
    for (it = obs.begin(); t <= obs.size() && it != obs.end(); ++it, ++t) {
        for (j = 0; j < m_nStates; ++j) {
            for (i = 0; i < m_nStates; ++i) {
                m_trelis[t][j] += m_trelis[t-1][i] * m_transitionProps[i][j];
            }
            m_trelis[t][j] *= m_outputProps[j][*it];
        }
    }
    return m_trelis;
}

Matrix<double>& LeftToRightHMM::backwardTrelis(Observation &obs)
{
    m_trelis = Matrix<double>(obs.size()+1, m_nStates, 0);
    // only state m_nStates-1 (= end state) has propability 1
    m_trelis[m_nStates-1][obs.size()] = 1;

    // beta_t(i) = sum{j=0}{N}{beta_{t+1}(j) * a_{i,j} * b_j(O_{t+1})}
    Observation::reverse_iterator it;
    size_t t = obs.size()-1, i, j;
    for (it = obs.rbegin(); t >= 0 && it != obs.rend(); ++it, --t) {
        for (i = 0; i < m_nStates; ++i) {
            for (j = 0; j < m_nStates; ++j) {
                m_trelis[t][i] += m_trelis[t+1][j] * m_transitionProps[i][j]
                                        * m_outputProps[j][*it];
            }
        }
    }
    return m_trelis;
}

---