from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
The state of the system is its position and velocity. The velocity is decreased by a damping factor. We want to infer the damping factor of the system.
# The state of the stystem
true_pos = 0.
true_vel = 1.
# Assume linear decrease of the speed of the system
true_damping = -1.5
We use simple euler integration to simulate the system.
# timestep
dt = 0.01
p = true_pos
v = true_vel
positions, velocities = [], []
while True:
# euler integration
p = p + v * dt
v = v + true_damping * dt
positions.append(p)
velocities.append(v)
if v <= 0:
break
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(positions, label="pos")
ax.plot(velocities, label="vel")
ax.legend()
ax.grid();
Now let's use PyMC to model this.
import pymc as pm
# let's treat the real positions as observations
observations = positions
# PRIORS
damping = pm.Normal("damping", mu=0, tau=1/20.**2)
#damping = pm.Uniform("damping", lower=-4, upper=4)
# we assume little system noise
tau_system_noise = (1 / 0.5**2)
# the state consist of pos and vel
# vel: we can't judge the initial velocity --> assume it's 0 with big std
vel_states = [pm.Normal("v0", mu=0, tau=(1 / 2**2))]
# pos: the first pos is just the observation
pos_states = [pm.Normal("s0", mu=observations[0], tau=tau_system_noise)]
for i in range(1, len(observations)):
new_vel = pm.Normal("v" + str(i),
mu=vel_states[-1] + damping * dt,
tau=(1/2.**2))
new_pos = pm.Normal("s" + str(i),
mu=pos_states[-1] + new_vel * dt,
tau=tau_system_noise)
vel_states.append(new_vel)
pos_states.append(new_pos)
# some observation noise
tau_observation_noise = (1 / 0.005**2)
obs = pm.Normal("obs", mu=pos_states, tau=tau_observation_noise, value=observations, observed=True)
mcmc = pm.MCMC([damping, vel_states, pos_states, obs])
mcmc.sample(10000, 5000)
pm.Matplot.plot(mcmc.get_node("damping"))
damping_samples = mcmc.trace("damping")[:]
print "damping -- mean:%f; std:%f" % (mean(damping_samples), std(damping_samples))
print "real damping -- %f" % true_damping
[-----------------100%-----------------] 10000 of 10000 complete in 181.4 secdamping -- mean:6.896805; std:16.894129 real damping -- -1.500000
v = np.vstack([mcmc.trace("v" + str(i))[:] for i in range(len(observations))])
s = np.vstack([mcmc.trace("s" + str(i))[:] for i in range(1,len(observations))])
plt.figure(figsize=(12,3))
plt.subplot(121)
plot(observations, lw = 2)
plot(s[:,::500]);
plt.subplot(122)
plot(v[:,::500]);