在PyMC中学习离散HMM参数

Question

我正在尝试使用PyMC学习一个简单的离散隐马尔可夫模型（HMM）的参数。我在模拟一个雨天和晴天的模型，这个模型可以在维基百科的HMM页面上找到。模型大概长这样：

我使用了以下的先验分布。

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的错误。

我不确定自己是否正确理解了这些回答。

Answer 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函数来获取合适的概率，然后把它作为分类的参数。

Answer 2

这里有一些想法：

• 晴天和雨天是互斥且完全的隐藏状态。为什么不把它们编码成一个天气变量，这个变量可以有两个值（一个代表晴天，一个代表雨天）呢？
• 在你的似然函数中，你似乎观察到了雨天/晴天。从你的模型图来看，这些似乎是隐藏变量，而不是观察到的变量（观察到的变量应该是“走路”、“购物”和“清洁”）。

在你的似然函数中，你需要对所有时间步骤 t 的观察值（走路、购物和清洁的值）进行求和，计算这些值在当前（采样的）雨天/晴天（即天气）下的对数概率。

学习

如果你想为模型学习参数，可能需要考虑切换到 PyMC3，这样可以更好地自动计算你的对数概率函数的梯度。不过在这种情况下（因为你选择了共轭先验），其实并不是特别必要。如果你不知道什么是共轭先验，或者需要了解一下，可以去维基百科查找共轭先验列表，那上面有很好的文章。

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

如果你对边际后验分布不感兴趣，而是想找到联合最大后验参数，你可以考虑使用期望最大化（EM）学习或模拟退火。这两种方法在 MCMC 框架内都应该能很好地工作。

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

• E（期望）步骤：运行 MCMC 链一段时间并收集样本
• M（最大化）步骤：将你的 Beta 和 Dirichlet 先验的超参数更新为这些样本就像它们是实际观察到的一样。如果你不知道怎么做，可以查一下 Beta 和 Dirichlet 分布。

我建议使用一个小的学习率，这样你就不会陷入第一个局部最优解（这时我们接近模拟退火）：与其把你从 MCMC 链生成的 N 个样本当作实际观察，不如把它们当作 K 个观察（K 的值远小于 N），通过将超参数的更新按 K/N 的学习率因子缩小。