1

我正在尝试使用 PyMC 学习简单离散 HMM 的参数。我正在对 HMM 上的Wiki页面中的雨晴模型进行建模。该模型如下所示：

我正在使用以下先验。

theta_start_state ~ beta(20,10)
theta_transition_rainy ~beta(8,2)
theta_transition_sunny ~beta(2,8)
theta_emission_rainy ~ Dirichlet(3,4,3)
theta_emission_sunny ~ Dirichlet(10,6,4)

最初，我使用此设置创建一个训练集，如下所示。

## Some not so informative priors!
# Prior on start state
theta_start_state = pm.Beta('theta_start_state',12,8)

# Prior on transition from rainy
theta_transition_rainy = pm.Beta('transition_rainy',8,2)

# Prior on transition from sunny
theta_transition_sunny = pm.Beta('transition_sunny',2,8)

# Prior on emission from rainy
theta_emission_rainy = pm.Dirichlet('emission_rainy',[3,4,3])

# Prior on emission from sunny
theta_emission_sunny = pm.Dirichlet('emission_sunny',[10,6,4])

# Start state
x_train_0 = pm.Categorical('x_0',[theta_start_state, 1-theta_start_state])

N = 100

# Create a train set for hidden states
x_train = np.empty(N, dtype=object)

# Creating a train set of observations
y_train = np.empty(N, dtype=object)

x_train[0] = x_train_0

for i in xrange(1, N):
    if x_train[i-1].value==0:
        x_train[i] = pm.Categorical('x_train_%d'%i,[theta_transition_rainy, 1- theta_transition_rainy])
    else:
        x_train[i] = pm.Categorical('x_train_%d'%i,[theta_transition_sunny, 1- theta_transition_sunny])


for i in xrange(0,N):
    if x_train[i].value == 0:
        # Rain
        y_train[i] = pm.Categorical('y_train_%d' %i, theta_emission_rainy)
    else:
        y_train[i] = pm.Categorical('y_train_%d' %i, theta_emission_sunny)

但是，我无法理解如何使用 PyMC 学习这些参数。我开始如下。

@pm.observed
def y(x=x_train, value =y_train):    
    N = len(x)    
    out = np.empty(N, dtype=object)
    for i in xrange(0,N):
        if x[i].value == 0:
            # Rain
            out[i] = pm.Categorical('y_%d' %i, theta_emission_rainy)
        else:
            out[i] = pm.Categorical('y_%d' %i, theta_emission_sunny)
    return out

可以在此处找到包含此代码的完整笔记本。

旁白：包含高斯 HMM 代码的要点真的很难理解！（未记录）

更新

根据以下答案，我尝试按如下方式更改我的代码：

@pm.stochastic(observed=True)
def y(value=y_train, hidden_states = x_train):
    def logp(value, hidden_states):
        logprob = 0 
        for i in xrange(0,len(hidden_states)):
            if hidden_states[i].value == 0:
                # Rain
                logprob = logprob + pm.categorical_like(value[i], theta_emission_rainy)
            else:
                # Sunny
                logprob = logprob + pm.categorical_like(value[i], theta_emission_sunny)
        return logprob

下一步将是创建一个模型，然后运行 ​​MCMC 算法。但是，上面编辑的代码也不起作用。它给出了一个ZeroProbability error.

我不确定我是否正确解释了答案。

4

2 回答 2

2

只是对此的一些想法：

  • Sunny 和 Rainy 是互斥且详尽的隐藏状态。为什么不将它们编码为单个分类天气变量，它可以采用两个值之一（编码为晴天、下雨天）？
  • 在您的似然函数中，您似乎观察到雨天/晴天。我在您的模型图中看到它的方式，这些似乎是隐藏的，而不是观察到的变量（即“步行”、“购物”和“清洁”）

在您的似然函数中，您需要在给定当前（采样）雨天/晴天（即天气）值的情况下，将（对于所有时间步长 t）观测值（分别为步行、购物和清洁）的对数概率求和同一时间步 t。

学习

如果您想学习模型的参数，您可能需要考虑切换到 PyMC3，它更适合自动计算 logp 函数的梯度。但是在这种情况下（因为您选择了共轭先验），这并不是必需的。如果您不知道共轭先验是什么，或者需要概述，请向维基百科询问共轭先验列表，它有一篇很棒的文章。

根据您想要做的事情，您有几个选择。如果你想从所有参数的后验分布中采样，只需像你一样指定你的 MCMC 模型，然后按下推理按钮，然后绘制并总结你感兴趣的参数的边际分布，你就完成了.

如果您对边际后验分布不感兴趣，而对寻找联合MAP 参数感兴趣，则可以考虑使用期望最大化 (EM) 学习或模拟退火。两者都应该在 MCMC 框架内工作得相当好。

对于 EM 学习，只需重复这些步骤直到收敛：

  • E（期望）步骤：运行MCMC链一段时间并收集样本
  • M（最大化）步骤：更新 Beta 和 Dirichlet Priors 的超参数，就好像这些样本是实际观察值一样。如果您不知道该怎么做，请查阅 Beta 和 Dirichlet Distributions。

我会使用一个小的学习率因子，这样您就不会陷入第一个局部最优值（现在我们正在接近模拟退火）：不要将您从 MCMC 链生成的 N 个样本视为实际观察值，而是将它们视为 K 个观察值（对于值 K << N）通过将超参数的更新按 K/N 的学习率因子缩小。

于 2014-05-06T21:05:42.750 回答
1

我首先想到的是您的可能性的返回值。PyMC 需要一个标量返回值，而不是列表/数组。您需要在返回之前对数组求和。

此外，当您使用 Dirichlet 作为分类的先验时，PyMC 会检测到这一点并填写最后一个概率。以下是我将如何编码您的x_train/y_train循环：

p = []
for i in xrange(1, N):
    # This will return the first variable if prev=0, and the second otherwise
    p.append(pm.Lambda('p_%i' % i, lambda prev=x_train[i-1]: (theta_transition_rainy, theta_transition_sunny)[bool(prev)]))
    x_train[i] = pm.Categorical('x_train_%i' % i, p[-1])

因此，您使用 Lambda 获取适当的概率，并将其用作分类的参数。

于 2014-05-06T03:36:24.187 回答