Posted By

fonnesbeck on 02/09/10


Tagged

states textmate python hidden PyMC Bayesian Markovian


Versions (?)

Who likes this?

1 person have marked this snippet as a favorite

adkatrit


Hidden Markov Model


 / Published in: Python
 

  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  

You need to login to post a comment.