Advanced Search

Hidden Markov Model

/ Published in: Python
Save to your folder(s)

Expand | Embed | Plain Text

Copy this code and paste it in your HTML 
  1. import numpy as np
  2. import pymc
  3. import pdb
  4.  
  5. def unconditionalProbability(Ptrans):
  6.    """Compute the unconditional probability for the states of a
  7.    Markov chain."""
  8.  
  9.    m = Ptrans.shape[0]
  10.  
  11.    P = np.column_stack((Ptrans, 1. - Ptrans.sum(axis=1)))
  12.  
  13.    I = np.eye(m)
  14.    U = np.ones((m, m))
  15.    u = np.ones(m)
  16.  
  17.    return np.linalg.solve((I - P + U).T, u)
  18.  
  19. data = np.loadtxt('test_data.txt',
  20.                  dtype=np.dtype([('state', np.uint8),
  21.                                  ('emission', np.float)]),
  22.                  delimiter=',',
  23.                  skiprows=1)
  24.  
  25. # Two state model for simplicity.
  26. N_states = 2
  27. N_chain = len(data)
  28.  
  29. # Transition probability stochastic
  30. theta = np.ones(N_states) + 1.
  31.  
  32. def Ptrans_logp(value, theta):
  33.    logp = 0.
  34.    for i in range(value.shape[0]):
  35.        logp = logp + pymc.dirichlet_like(value[i], theta)
  36.    return logp
  37.  
  38. def Ptrans_random(theta):
  39.    return pymc.rdirichlet(theta, size=len(theta))
  40.  
  41. Ptrans = pymc.Stochastic(logp=Ptrans_logp,
  42.                         doc='Transition matrix',
  43.                         name='Ptrans',
  44.                         parents={'theta': theta},
  45.                         random=Ptrans_random)
  46.  
  47. #Hidden states stochastic
  48. def states_logp(value, Ptrans=Ptrans):
  49.  
  50.     if sum(value>1):
  51.         return -np.inf
  52.  
  53.     P = np.column_stack((Ptrans, 1. - Ptrans.sum(axis=1)))
  54.  
  55.     Pinit = unconditionalProbability(Ptrans)
  56.  
  57.     logp = pymc.categorical_like(value[0], Pinit)
  58.  
  59.     for i in range(1, len(value)):
  60.        try:
  61.           logp = logp + pymc.categorical_like(value[i], P[value[i-1]])
  62.        except:
  63.           pdb.set_trace()
  64.  
  65.     return logp
  66.  
  67. def states_random(Ptrans=Ptrans, N_chain=N_chain):
  68.    P = np.column_stack((Ptrans, 1. - Ptrans.sum(axis=1)))
  69.  
  70.    Pinit = unconditionalProbability(Ptrans)
  71.  
  72.    states = np.empty(N_chain, dtype=np.uint8)
  73.  
  74.    states[0] = pymc.rcategorical(Pinit)
  75.  
  76.    for i in range(1, N_chain):
  77.        states[i] = pymc.rcategorical(P[states[i-1]])
  78.  
  79.    return states
  80.  
  81. states = pymc.Stochastic(logp=states_logp,
  82.                         doc='Hidden states',
  83.                         name='states',
  84.                         parents={'Ptrans': Ptrans},
  85.                         random=states_random,
  86.                         dtype=np.uint8)
  87.  
  88. # Gaussian emission parameters
  89. mu = pymc.Normal('mu', 0., 1.e-6, value=np.random.randn(N_states))
  90. sigma = pymc.Uniform('sigma', 0., 100.,
  91. value=np.random.rand(N_states))
  92.  
  93. y = pymc.Normal('y', mu[states], 1./sigma[states]**2,
  94. value=data['emission'], observed=True)

Report this snippet


Comments

RSS Icon Subscribe to comments

You need to login to post a comment.