Viterbi Algorithm implemented incorrectly

[mo-kli]

@DonatoMeoli @antmarakis I think the backtracking in the Viterbi algorithm added in #1099 is implemented incorrectly. It becomes obvious when seeing that the way it is implemented, the backtracking could also be applied as a forward algorithm:

    for i in range(t, -1, -1):
        path[i - 1] = max(m[i - 1])

When maximizing over the previous m

m[i] = element_wise_product(HMM.sensor_dist(ev[i + 1]),
                 [max(element_wise_product(HMM.transition_model[0], m[i - 1])),
                  max(element_wise_product(HMM.transition_model[1], m[i - 1]))])

there should be kept track of the maximizing element, e.g.

backtracking_graph.append([argmax(element_wise_product(HMM.transition_model[0], m[i - 1])),
                           argmax(element_wise_product(HMM.transition_model[1], m[i - 1]))])`

Afterwards, backtracking can be implemented as

i_max = argmax(m[-1])
    for i in range(t-1, -1, -1):
        ml_sequence_2[i-1] = m[i-1][i_max]
        if i > 0:
            i_max = backtracking_graph[i-1][i_max]

In my opinion, the complete function should look like:

from numpy import argmax

def viterbi(HMM, ev):
    """[Equation 15.11]
    Viterbi algorithm to find the most likely sequence. Computes the best path,
    given an HMM model and a sequence of observations."""
    t = len(ev)
    ev.insert(0, None)

    m = [[0.0, 0.0] for _ in range(len(ev) - 1)]

    # the recursion is initialized with m1 = forward(P(X0), e1)
    m[0] = forward(HMM, HMM.prior, ev[1])
    backtracking_graph = []

    for i in range(1, t):
        m[i] = element_wise_product(HMM.sensor_dist(ev[i + 1]),
                                    [max(element_wise_product(HMM.transition_model[0], m[i - 1])),
                                     max(element_wise_product(HMM.transition_model[1], m[i - 1]))])
        backtracking_graph.append([argmax(element_wise_product(HMM.transition_model[0], m[i - 1])),
                                   argmax(element_wise_product(HMM.transition_model[1], m[i - 1]))])

    ml_sequence = [0.0] * (len(ev) - 1)

    # the construction of the most likely sequence starts in the final state with the largest probability,
    # and runs backwards; the algorithm needs to store for each xt its maximizing predecessor xt-1
    i_max = argmax(m[-1])
    for i in range(t-1, -1, -1):
        ml_sequence[i] = m[i][i_max]
        if i > 0:
            i_max = backtracking_graph[i-1][i_max]

    return ml_sequence

To see an example how backtracking is applied, see also here,

[dmeoli]

Ok @Tsovet it seems good. Could you proceed with a pull request?

[mo-kli]

@DonatoMeoli I just created a pull request

[dmeoli]

Well done! Waitin' @antmarakis for the merge