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

damping的值由先验决定。即使我将之前的改变为制服或其他什么，它仍然是这样的。

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

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

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

EDIT2：@Chris answer没有帮助。我不能使用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())

Answer 2

假设您知道模型的结构-您正在进行参数估计，而不是系统识别-您可以将PyMC模型构建为回归模型，将未知阻尼、初始位置和初始速度作为参数以及位置数组和观测值。

也就是说，PM类表示点质量系统：

pm = PM(true_damping)

positions, velocities = pm.integrate(true_pos, true_vel, N, dt)

# Assume little system noise
std_system_noise = 0.05
tau_system_noise = 1.0/std_system_noise**2

# Treat the real positions as observations
observations = positions + np.random.randn(N,)*std_system_noise

# Damping is modelled with a Uniform prior
damping = mc.Uniform("damping", lower=-4.0, upper=4.0, value=-0.5)

# Initial position & velocity unknown -> assume Uniform priors
init_pos = mc.Uniform("init_pos", lower=-1.0, upper=1.0, value=0.5)
init_vel = mc.Uniform("init_vel", lower=0.0, upper=2.0, value=1.5)

@mc.deterministic
def det_pos(d=damping, pi=init_pos, vi=init_vel):
    # Apply damping, init_pos and init_vel estimates and integrate
    pm.damping = d.item()
    pos, vel = pm.integrate(pi, vi, N, dt)
    return pos

# Standard deviation is modelled with a Uniform prior
std_pos = mc.Uniform("std", lower=0.0, upper=1.0, value=0.5)

@mc.deterministic
def det_prec_pos(s=std_pos):
    # Precision, based on standard deviation
    return 1.0/s**2

# The observations are based on the estimated positions and precision
obs_pos = mc.Normal("obs", mu=det_pos, tau=det_prec_pos, value=observations, observed=True)

# Create the model and sample
model = mc.Model([damping, init_pos, init_vel, det_prec_pos, obs_pos])
mcmc = mc.MCMC(model)
mcmc.sample(50000, 25000)

完整的清单在这里：https://gist.github.com/stuckeyr/7762371

增加N和减少dt将显著改善您的估计。

Answer 3

你说的不合理是什么意思？它们是不是向前一个方向缩小？阻尼似乎有一个非常紧密的变化--如果你给它一个更分散的先验呢？

另外，您可以尝试在*_states数组上使用AdaptiveMetropolis采样器：

my_model.use_step_method(AdaptiveMetropolis, my_model.vel_states)

对于相关变量，它有时会混合得更好，就像这些变量可能是的那样。