python - PyMC:分层隐马尔可夫模型

Question

这是 PyMC: Parameter estimation in a Markov system 的后续内容

我有一个由每个时间步的位置和速度定义的系统。系统的行为定义为:

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

阻尼的值由先验决定。即使我将之前的值更改为 Uniform 或其他什么，情况仍然如此。

我做错了什么？它与前面的示例非常相似，只是多了一层。

此问题的完整 IPython 笔记本可在此处找到:http://nbviewer.ipython.org/github/sotte/random_stuff/blob/master/PyMC%20-%20HMM%20Dynamic%20System.ipynb

[编辑:一些说明和采样代码。]

[EDIT2:@Chris 的回答没有帮助。我无法使用 AdaptiveMetropolis 因为 *_states 似乎不是模型的一部分。]

Answer 1

该模型存在一些问题，请再次查看。首先也是最重要的，您没有将所有 PyMC 对象添加到模型中。您仅添加了[damping, obs]。您应该将所有 PyMC 节点传递给模型。

另请注意，您无需同时调用 Model 和 MCMC。这很好:

model = pm.MCMC([damping, obs, vel_states, pos_states])

PyMC 的最佳工作流程是将模型保存在与运行逻辑分开的文件中。这样，您只需导入模型并将其传递给 MCMC:

import my_model

model = pm.MCMC(my_model)

或者，您可以将模型编写为函数，返回locals(或vars)，然后调用该函数作为MCMC的参数>。例如:

def generate_model():

    # put your model definition here

    return locals()

model = pm.MCMC(generate_model())