1

在阅读了 Cam Davidson-Pilon 的Probabilistic Programming & Bayesian Methods for Hackers之后，我决定尝试使用 PyMC解决隐马尔可夫模型 (HMM)学习问题。到目前为止，代码不配合，但是通过排查，感觉已经缩小了问题的根源。

将代码分解成更小的块并专注于 t=0 时的初始概率和发射概率，我能够学习 t=0 时单个状态的发射/观察参数。但是，一旦我添加另一个状态（总共两个状态），无论数据输入如何，参数学习的结果都是相同的（并且不正确）。所以，我觉得我一定在代码部分做错了@pm.deterministic，这不允许我从Init初始概率函数中采样。

通过这部分代码，我的目标是学习初始概率 p_bern和发射概率 p_0，p_1分别对应于状态 0 和 1。发射取决于状态，这是我试图用我的@pm.deterministic函数来表达的。我可以在这个确定性函数中使用“if”语句吗？这似乎是问题的根源。

# This code is to test the ability to discern between two states with emissions

import numpy as np
import pymc as pm
from matplotlib import pyplot as plt

N = 1000
state = np.zeros(N)
data = np.zeros(shape=N)

# Generate data
for i in range(N):
    state[i] = pm.rbernoulli(p=0.3)
for i in range(N):
    if state[i]==0:
        data[i] = pm.rbernoulli(p=0.4)
    elif state[i]==1:
        data[i] = pm.rbernoulli(p=0.8)

# Prior on probabilities
p_bern = pm.Uniform("p_S", 0., 1.)
p_0 = pm.Uniform("p_0", 0., 1.)
p_1 = pm.Uniform("p_1", 0., 1.)

Init = pm.Bernoulli("Init", p=p_bern) # Bernoulli node

@pm.deterministic
def p_T(Init=Init, p_0=p_0, p_1=p_1, p_bern=p_bern):
    if Init==0:
        return p_0
    elif Init==1:
        return p_1

obs = pm.Bernoulli("obs", p=p_T, value=data, observed=True)
model = pm.Model([obs, p_bern, p_0, p_1])
mcmc = pm.MCMC(model)
mcmc.sample(20000, 10000)
pm.Matplot.plot(mcmc)

我已经尝试了以下无济于事：

  1. 使用@pm.potential装饰器创建联合分布
  2. 更改我的位置的Init位置（您可以在代码中看到我的注释，我不确定将它放在哪里）
  3. 使用@pm.stochastic与此类似的

编辑：根据克里斯的建议，我已经将伯努利节点移到了确定性之外。我还将代码更新为更简单的模型（伯努利观察而不是多项式），以便于排除故障。

感谢您的时间和关注。任何反馈都会受到热烈欢迎。另外，如果我缺少任何信息，请告诉我！

4

2 回答 2

2

根据您提供的更新信息，以下是一些有效的代码：

import numpy as np
import pymc as pm
from matplotlib import pyplot as plt

N = 1000
state = np.zeros(N)
data = np.zeros(shape=N)

# Generate data
state = pm.rbernoulli(p=0.3, size=N)
data = [int(pm.rbernoulli(0.8*s or 0.4)) for s in state]

# Prior on probabilities
p_S = pm.Uniform("p_S", 0., 1.)
p_0 = pm.Uniform("p_0", 0., 1.)
p_1 = pm.Uniform("p_1", 0., 1.)

# Use values of Init as indices to probabilities
Init = pm.Bernoulli("Init", p=p_S, size=N) # Bernoulli node

p_T = pm.Lambda('p_T', lambda p_0=p_0, p_1=p_1, i=Init: np.array([p_0, p_1])[i.astype(int)])

obs = pm.Bernoulli("obs", p=p_T, value=data, observed=True)
model = pm.MCMC(locals())
model.sample(20000, 10000)
model.summary()

请注意，在数据生成步骤中，我使用状态来索引适当的真实概率。我基本上在p_T. 它似乎工作得相当好，但请注意，根据初始化的位置， 和 的两个值p_0最终p_1可能对应于任何一个真实值（没有什么可以限制一个大于另一个）。因此， 的值p_S最终可以作为真实状态概率的补数。

于 2015-08-30T21:36:25.713 回答
2

我会将这种随机性从确定性中移出。确定性节点的值应该完全由其父节点的值决定。将随机变量埋在节点中违反了这一点。

为什么点创建一个伯努利节点，并将其作为参数传递给确定性？

于 2015-08-24T01:22:35.050 回答