4

A Simple Markov Chain

Let's say we want to estimate parameters of a system such that we can predict the state of the system at time step 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.

bocebir
  • 13
  • 1
  • 4
Stefan
  • 1,246
  • 1
  • 9
  • 13

1 Answers1

3

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.

Chris Fonnesbeck
  • 4,143
  • 4
  • 29
  • 30
  • Thanks, it works. Using list comprehensions does not work though, because you can't easily access the last added element in a list comprehension. I've updated my toy problem and incorporated your solution: http://nbviewer.ipython.org/github/sotte/random_stuff/blob/master/PyMC%20-%20Simple%20Markov%20Chain.ipynb – Stefan Nov 28 '13 at 11:06