# Hidden Markov Model

Apr 29th, 2012
98
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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.
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)
RAW Paste Data