Calculating HMM observation probability

Question

I'm putting this one in stackoverflow rather than math.stackexchange, seeing as I'm trying a programming approach rather than a math approach.

I have 3 matrices, a transition trans, an emission emiss (or observation error) and a state state.

To go with this, I also have a series of observations obs

My approach is rather simple:

prob = 1                       # Probability of sequence is 1 to start with
for o in obs:                  # For each observation:
    p_mult = dot(emiss, state) # Get the probability of each observation,
    p_mult = p_mult.get_elem(o)# select the corresponding observation
    prob = prob * p_mult       # and multiply that with the total probability
    state = dot(trans, state)  # Last, change the state using the transition matrix
print(prob)                    # Print answer

Where dot(x,y) is the dot product of two matrices (eg {1x4}*{4x3}->{1x3}) and x.get_elem(y) takes the yth element of the vector x .

For some reason this does not seem to work, as the probabilities I'm calculating is not matching that of others. Can someone give me a hint of what is wrong with this reasoning?