from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
%pylab inline
Couldn't import dot_parser, loading of dot files will not be possible. Populating the interactive namespace from numpy and matplotlib
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 $t$. PyMC should be able to deal with this easily.
Let our system be 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_pos, true_vel = 0, .7
true_positions = [true_vel * step for step in range(100)]
We assume that we have sligh 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$.
import random
# we're using `some_tau` for the noise throughout the example.
# this should be replaced with something more meaningful.
some_tau = 1 / .5**2
# PRIORS
# we don't know too much about the velocity, might be pos. or neg.
vel = pm.Normal("vel", mu=0, tau=some_tau)
# MODEL
# next_state = prev_state + vel (and some gaussian noise)
# That means that each state depends on the prev_state and the vel.
# We save the states in a list.
states = [pm.Normal("s0", mu=true_positions[0], tau=some_tau)]
for i in range(1, len(true_positions)):
states.append(pm.Normal(name="s" + str(i),
mu=states[-1] + vel,
tau=some_tau))
# observation with gaussian noise
obs = pm.Normal("obs", mu=states, tau=some_tau, value=true_positions, observed=True)
model = pm.Model([vel, obs])
mcmc = pm.MCMC(model)
mcmc.sample(10000, 5000)
[****************100%******************] 10000 of 10000 complete
vel_samples = mcmc.trace("vel")[:]
print "vel -- mean:%f; std:%f" % (mean(vel_samples), std(vel_samples))
pm.Matplot.plot(mcmc.get_node("vel"))
vel -- mean:0.694063; std:0.051844 Plotting vel