我有一个系统,它由每个时间步的位置和速度定义。系统的行为定义为:
vel = vel + damping * dt
pos = pos + vel * dt
所以,这是我的 PyMC 模型。估计vel
,pos
也是最重要的damping
。
# PRIORS
damping = pm.Normal("damping", mu=-4, tau=(1 / .5**2))
# we assume some system noise
tau_system_noise = (1 / 0.1**2)
# the state consist of (pos, vel); save in lists
# vel: we can't judge the initial velocity --> assume it's 0 with big std
vel_states = [pm.Normal("v0", mu=-4, tau=(1 / 2**2))]
# pos: the first pos is just the observation
pos_states = [pm.Normal("p0", 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=tau_system_noise)
vel_states.append(new_vel)
pos_states.append(
pm.Normal("s" + str(i),
mu=pos_states[-1] + new_vel * dt,
tau=tau_system_noise)
)
# we assume some observation noise
tau_observation_noise = (1 / 0.5**2)
obs = pm.Normal("obs", mu=pos_states, tau=tau_observation_noise, value=observations, observed=True)
这就是我运行采样的方式:
mcmc = pm.MCMC([damping, obs, vel_states, pos_states])
mcmc.sample(50000, 25000)
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
的值damping
由先验支配。即使我将之前的内容更改为 Uniform 或其他内容,情况仍然如此。
我究竟做错了什么?它与前面的示例非常相似,只是多了一层。
此问题的完整 IPython 笔记本可在此处获得:http: //nbviewer.ipython.org/github/sotte/random_stuff/blob/master/PyMC%20-%20HMM%20Dynamic%20System.ipynb
[编辑:一些说明和采样代码。]
[EDIT2:@Chris 的回答没有帮助。我无法使用AdaptiveMetropolis
,因为 *_states 似乎不是模型的一部分。]