A Simple Markow Chain

Let's say we want to estimate parameters of a system such that we can predict the state of the system at timestep t+1 given the state at timestep t. PyMC should be able to deal with this easily.

Let our toy system consist of a moving object in a 1D world. The state is the position of the object. We want to estimate the latent variable/the speed of the object. The next state depends on the previous state and the latent variable the speed.

# define the system and the data
true_vel = .2
true_pos = 0
true_positions = [.2 * step for step in range(100)]

We assume that we have some noise in our observation (but that does not matter here).

The question is: how do I model the dependency of the next state on the current state. I could supply the transition function a parameter idx to access the position at time t and then predict the position at time t+1.

vel = pymc.Normal("pos", 0, 1/(.5**2))
idx = pymc.DiscreteUniform("idx", 0, 100, value=range(100), observed=True)

@pm.deterministic
def transition(positions=true_positions, vel=vel, idx=idx):
    return positions[idx] + vel

# observation with gaussian noise
obs = pymc.Normal("obs", mu=transition, tau=1/(.5**2))

However, the index seems to be an array which is not suitable for indexing. There is probably a better way to access the previous state.

The easiest way is to generate a list, and allow PyMC to deal with it as a Container. There is a relevant example on the PyMC wiki. Here is the relevant snippet:

# Lognormal distribution of P's
Pmean0 = 0.
P_0 = Lognormal('P_0', mu=Pmean0, tau=isigma2, trace=False, value=P_inits[0])
P = [P_0]

# Recursive step
for i in range(1,nyears):
    Pmean = Lambda("Pmean", lambda P=P[i-1], k=k, r=r: log(max(P+r*P*(1-P)-k*catch[i-1],0.01)))
    Pi = Lognormal('P_%i'%i, mu=Pmean, tau=isigma2, value=P_inits[i], trace=False)
    P.append(Pi)

Notice how the mean of the current Lognormal is a function of the last one? Not elegant, using list.append and all, but you can use a list comprehension instead.