可以将“if”语句用于PyMC确定性函数吗？

在阅读Cam Davidson-Pilon的Probabilistic Programming & Bayesian Methods for Hackers后，我决定尝试用隐马尔可夫模型（HMM）与PyMC学习问题。到目前为止，代码并不合作，但通过故障排除，我觉得我已经缩小了问题的根源。可以将“if”语句用于PyMC确定性函数吗？

将代码分解成更小的块并关注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)

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

使用@pm.potential装饰以create a joint distribution
更改我的Init位置的位置（你可以看到我的评论的代码，我我不确定放在哪里）
使用@pm.stochastic类似于this

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

谢谢你的时间和关注。任何反馈都会受到热烈的欢迎。另外，如果我错过了任何信息，请让我知道！

来源

2015-08-14 Richard

我会将这种随机性排除在确定性之外。确定性节点的价值应该完全由其父母的价值决定。掩埋在节点中的随机变量违反了这一点。

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

来源

2015-08-24 01:22:35

嗨@Chris，非常感谢您看看我的代码。你指的是将'Init = pm.rbernoulli（p = p_bern）'移出确定性，并将其作为参数传递，如'Init == 0：'？如果是这样，我已经尝试了这种方法，正如你可以在前面看到的'ed'行中看到的那样，得到相同的结果。或者你的意思是别的吗？ – Richard

就是这条线。它导致确定性节点不再是确定性的。 –

因此在逻辑上这意味着在确定性中使用if语句没有任何错误？ – Richard

根据您提供的最新信息，这里是一些代码，工作原理：在数据生成步骤

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-30 21:36:25

嗨克里斯，谢谢你的帮助。您使用状态进行索引的方法非常有用。实际上，您的数据生成和我的有什么不同？或者这是效率问题？另外，我发现你使用数组的方法来指定'p_T'非常酷。你为什么选择这个方法而不是使用'if'语句？ – Richard

可以将“if”语句用于PyMC确定性函数吗？

回答

相关问题