/ Published in: Python
Expand |
Embed | Plain Text
Copy this code and paste it in your HTML
import numpy as np import pymc import pdb def unconditionalProbability(Ptrans): """Compute the unconditional probability for the states of a Markov chain.""" m = Ptrans.shape[0] P = np.column_stack((Ptrans, 1. - Ptrans.sum(axis=1))) I = np.eye(m) U = np.ones((m, m)) u = np.ones(m) return np.linalg.solve((I - P + U).T, u) data = np.loadtxt('test_data.txt', dtype=np.dtype([('state', np.uint8), ('emission', np.float)]), delimiter=',', skiprows=1) # Two state model for simplicity. N_states = 2 N_chain = len(data) # Transition probability stochastic theta = np.ones(N_states) + 1. def Ptrans_logp(value, theta): logp = 0. for i in range(value.shape[0]): logp = logp + pymc.dirichlet_like(value[i], theta) return logp def Ptrans_random(theta): return pymc.rdirichlet(theta, size=len(theta)) Ptrans = pymc.Stochastic(logp=Ptrans_logp, doc='Transition matrix', name='Ptrans', parents={'theta': theta}, random=Ptrans_random) #Hidden states stochastic def states_logp(value, Ptrans=Ptrans): if sum(value>1): return -np.inf P = np.column_stack((Ptrans, 1. - Ptrans.sum(axis=1))) Pinit = unconditionalProbability(Ptrans) logp = pymc.categorical_like(value[0], Pinit) for i in range(1, len(value)): try: logp = logp + pymc.categorical_like(value[i], P[value[i-1]]) except: pdb.set_trace() return logp def states_random(Ptrans=Ptrans, N_chain=N_chain): P = np.column_stack((Ptrans, 1. - Ptrans.sum(axis=1))) Pinit = unconditionalProbability(Ptrans) states = np.empty(N_chain, dtype=np.uint8) states[0] = pymc.rcategorical(Pinit) for i in range(1, N_chain): states[i] = pymc.rcategorical(P[states[i-1]]) return states states = pymc.Stochastic(logp=states_logp, doc='Hidden states', name='states', parents={'Ptrans': Ptrans}, random=states_random, dtype=np.uint8) # Gaussian emission parameters mu = pymc.Normal('mu', 0., 1.e-6, value=np.random.randn(N_states)) sigma = pymc.Uniform('sigma', 0., 100., value=np.random.rand(N_states)) y = pymc.Normal('y', mu[states], 1./sigma[states]**2, value=data['emission'], observed=True)