1

这是PyMC 的后续：马尔可夫系统中的参数估计

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

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 似乎不是模型的一部分。]

4

4 回答 4

1

假设您知道模型的结构——您正在进行参数估计，而不是系统识别——您可以将 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 将显着改善您的估计。

于 2013-12-03T01:46:03.413 回答
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())
于 2013-12-02T15:15:29.260 回答
0

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

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

my_model.use_step_method(AdaptiveMetropolis, my_model.vel_states)

它有时会更好地混合相关变量，因为这些可能是。

于 2013-11-30T20:52:09.230 回答
0

我认为您的初始方法很好并且应该可以工作，除了“obs”变量未包含在提供给 MCMC 的节点列表中（请参阅笔记本中的 In[10]）。包含此变量后，MCMC 求解器运行良好，并强制执行模型指定的条件依赖项。我想重复一下 Chris 的观点，即最好在不同的文件或函数下定义模型以避免此类错误。

您没有得到正确结果的原因是您的先验是任意选择的，并且在某些情况下，这些值使得模型很难正确混合以收敛。你的玩具问题试图估计一个阻尼值，使位置收敛到观察位置的向量。为此，您的模型应该能够灵活地在较宽的范围内选择速度和阻尼值，以便在从一个时间步长到下一个时间步长时可以纠正位置/速度的随机误差。否则，由于您的欧拉积分方案，错误只会不断传播。我认为 Chris 在建议选择更分散的先验时提到了同样的事情。

我建议使用每个 Normal 变量的 tau 值。例如，我更改了以下值：

damping = pm.Normal("damping", mu=0, tau=1/20.**2) # was tau=1/2.**2

new_vel = pm.Normal("v" + str(i),
                    mu=vel_states[-1] + damping * dt,
                    tau=(1/2.**2)) # was tau=tau_system_noise=(1 / 0.5**2)

tau_observation_noise = (1 / 0.005**2) # was 1 / 0.5**2

您可以在此处查看修改后的文件。

底部的图表显示，这些位置确实在趋同。速度无处不在。阻尼的估计平均值为 6.9，与 -1.5 有很大不同。也许您可以通过为先验选择适当的值来获得更好的估计。

于 2013-12-04T01:00:53.797 回答