PyMC: parameter estimation in the Markov system

Question

Simple Markow Chain

Say we want to evaluate the parameters of the system so that we can predict the state of the system in time t + 1 under the condition of state in time t. PyMC should be able to handle this easily.

Let our toy system consist of a moving object in a 1D world. State is the position of the object. We want to evaluate the hidden variable / object speed. The next state depends on the previous state and the hidden variable 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 it does not matter here).

The question arises: how can I model the dependence of the next state on the current state. I could provide the transition function with the idx parameter 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 is represented as an array, which is not suitable for indexing. There is probably a better way to access the previous state.

Answer 1

The easiest way is to create a list and let PyMC treat it as a Container. There is a corresponding 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 average of the current Lognormal is a function of the latter? Not elegant using list.append and that’s it, but you can use list comprehension instead.