

In [3]:

    
observed_states = np.tile(np.array([0,30,50,100,100,75,200,20,10,0]),1000)
mus = np.array([0,100,50])
sigmas = np.array([20,20,20])
trans_prob = np.array([[.9,.05,.05],[.1,.8,.1],[.2,.2,.6]])
start_prob = -np.log(np.array([.3,.3,.3]))

def norm_pdf(x,params):
    mu,sigma = params
    return -np.log(np.exp(-(x-mu)**2/2/(sigma**2))/np.sqrt(2*np.pi)/sigma)

def state_probabilities(observed_states,pdf,params):
    state_probs = []
    for param in params:
        state_probs.append(pdf(observed_states,param))
    return np.array(state_probs).T

def get_states(observed_states,pdf,pdf_params,start_probs,trans_probs):
    state_probs = state_probabilities(observed_states,pdf,pdf_params)
    n_obs = state_probs.shape[0]
    n_states = state_probs.shape[1]
    start_probs = -np.log(start_probs)
    trans_probs = -np.log(trans_probs)
    
    init = start_probs + state_probs[0,:]
    prevs = np.zeros(state_probs.shape,dtype=np.int)
    probs = np.zeros(state_probs.shape)
    values = state_probs[1:,:,np.newaxis] + trans_probs[np.newaxis,:,:]
    probs[0] = init
    prevs[0] = np.nan
    
    for t,options in enumerate(values):
        possibilities = probs[t] + options
        probs[t+1] = np.min(possibilities,axis=1)
        prevs[t+1] = np.argmin(possibilities,axis=1)
    
    best = [np.argmin(probs[-1])]
    for prevs in prevs[-1:0:-1]:
        best.append(prevs[best[-1]])
    return np.array(best[::-1])
    

states = get_states(observed_states,norm_pdf,zip(mus,sigmas),start_prob,trans_prob)
print states









    



[0 2 2 ..., 0 0 0]