可以在 PyMC 确定性函数中使用“if”语句吗？

Question

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

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

通过这部分代码，我的目标是学习对应于状态 0 和 1 的 初始概率 p_bern 和 发射概率 p_0 和 p_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. 对create a joint distribution 使用@pm.potential 装饰器
2. 更改我的 Init 位置的位置（您可以在代码中看到我的注释，我不确定该放在哪里）
3. 使用类似于this 的@pm.stochastic

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

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

Answer 1

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

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

Answer 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 的值最终可以作为真实状态概率的补码。