I am trying to figure out how to properly make a discrete state Markov chain model with pymc
.
As an example (view in nbviewer), lets make a chain of length T=10 where the Markov state is binary, the initial state distribution is [0.2, 0.8] and that the probability of switching states in state 1 is 0.01 while in state 2 it is 0.5
import numpy as np
import pymc as pm
T = 10
prior0 = [0.2, 0.8]
transMat = [[0.99, 0.01], [0.5, 0.5]]
To make the model, I make an array of state variables and an array of transition probabilities that depend on the state variables (using the pymc.Index function)
states = np.empty(T, dtype=object)
states[0] = pm.Categorical('state_0', prior0)
transPs = np.empty(T, dtype=object)
transPs[0] = pm.Index('trans_0', transMat, states[0])
for i in range(1, T):
states[i] = pm.Categorical('state_%i' % i, transPs[i-1])
transPs[i] = pm.Index('trans_%i' %i, transMat, states[i])
Sampling the model shows that the states marginals are what they should be (compared with model built with Kevin Murphy's BNT package in Matlab)
model = pm.MCMC([states, transPs])
model.sample(10000, 5000)
[np.mean(model.trace('state_%i' %i)[:]) for i in range(T)]
prints out:
[-----------------100%-----------------] 10000 of 10000 complete in 7.5 sec
[0.80020000000000002,
0.39839999999999998,
0.20319999999999999,
0.1118,
0.064199999999999993,
0.044600000000000001,
0.033000000000000002,
0.026200000000000001,
0.024199999999999999,
0.023800000000000002]
My question is - this does not seem like the most elegant way to build a Markov chain with pymc. Is there a cleaner way that does not require the array of deterministic functions?
My goal is to write a pymc based package for more general dynamic Bayesian networks.