我可以使用什么 pymc3 Monte-Carlo 步进器进行自定义分类分布？

Question

我正在致力于在 pymc3 中实现隐藏马尔可夫链。我在实现隐藏状态方面已经取得了很大进展。下面，我展示了一个简单的二态马尔可夫链：

import numpy as np
import pymc3 as pm
import theano.tensor as tt

# Markov chain sample with 2 states that was created
# to have prob 0->1 = 0.1 and prob 1->0 = 0.3
sample = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0,
   0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
   1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,
   0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0],
dtype=np.uint8)

我现在正在定义描述状态的 class。作为输入，我需要知道 P1 从状态 0 移动到状态 1 的概率，以及 P2 从 1->0 移动的概率。我还需要知道第一个状态为 0 的概率 PA。

class HMMStates(pm.Discrete):
    """
    Hidden Markov Model States
    Parameters
    ----------
    P1 : tensor
        probability to remain in state 1
    P2 : tensor
        probability to move from state 2 to state 1

    """

    def __init__(self, PA=None, P1=None, P2=None,
                 *args, **kwargs):
        super(HMMStates, self).__init__(*args, **kwargs)
        self.PA = PA
        self.P1 = P1
        self.P2 = P2
        self.mean = 0.
        self.mode = tt.cast(0,dtype='int64')

    def logp(self, x):
        PA = self.PA
        P1 = self.P1
        P2 = self.P2

        # now we need to create an array with probabilities
        # so that for x=A: PA=P1, PB=(1-P1)
        # and for x=B: PA=P2, PB=(1-P2)

        choice = tt.stack((P1,P2))
        P = choice[x[:-1]]

        x_i = x[1:]

        ou_like = pm.Categorical.dist(P).logp(x_i)
        return pm.Categorical.dist(PA).logp(x[0]) + tt.sum(ou_like)

我对我在 theano google 组学到的高级索引忍者技巧感到非常自豪。您也可以使用 tt.switch 实现相同的功能。我不太确定的是 self.mode。我只是给它 0 以避免测试值错误。下面是如何在模型中使用 class 来测试它是否有效。在这种情况下，状态不是隐藏的，而是观察到的。

with pm.Model() as model:
    # 2 state model
    # P1 is probablility to stay in state 1
    # P2 is probability to move from state 2 to state 1
    P1 = pm.Dirichlet('P1', a=np.ones(2))
    P2 = pm.Dirichlet('P2', a=np.ones(2))

    PA = pm.Deterministic('PA',P2/(P2+1-P1))

    states = HMMStates('states',PA,P1,P2, observed=sample)

    start = pm.find_MAP()

    trace = pm.sample(5000, start=start)

输出很好地再现了数据。在下一个模型中，我将展示这个问题。在这里，我不直接观察状态，而是观察添加了一些高斯噪声的状态（因此是隐藏状态）。如果您运行使用 Metropolis 步进器的模型，那么它会因索引错误而崩溃，我将其追溯到与在分类分布上使用 Metropolis 步进器相关的问题。不幸的是，唯一适用于我的 class 的步进器是 CategoricalGibbsMetropolis 步进器，但它拒绝与我的 class 一起使用，因为它不是明确的分类分布。

gauss_sample = sample*1.0 + 0.1*np.random.randn(len(sample))
from scipy import optimize
with pm.Model() as model2:
    # 2 state model
    # P1 is probablility to stay in state 1
    # P2 is probability to move from state 2 to state 1
    P1 = pm.Dirichlet('P1', a=np.ones(2))
    P2 = pm.Dirichlet('P2', a=np.ones(2))

    S = pm.InverseGamma('S',alpha=2.1, beta=1.1)

    PA = pm.Deterministic('PA',P2/(P2+1-P1))

    states = HMMStates('states',PA,P1,P2, shape=len(gauss_sample))

    emission = pm.Normal('emission',
                         mu=tt.cast(states,dtype='float64'),
                         sd=S,
                         observed = gauss_sample)

    start2 = pm.find_MAP(fmin=optimize.fmin_powell)
    step1 = pm.Metropolis(vars=[P1, P2, S, PA, emission])
    step2 = pm.ElemwiseCategorical(vars=[states], values=[0,1])
    trace2 = pm.sample(10000, start=start, step=[step1,step2])

ElemwiseCategorical 使其成为运行，但没有为我的状态分配正确的值。状态要么全为0，要么全为1。

如何告诉 ElemwiseCategorial 分配 1 和 0 的状态向量，或者如何让 CategorialGibbsMetropolis 将我的分布识别为分类分布。这一定是自定义发行版的常见问题。

Answer 1

由于没有人回答我的问题，所以我自己回答了。我使用的技巧是 Chris Fonnesbeck 在我打开问题的 pymc3 github 上建议的。他建议subclass pm.Categorical.

class HMMStates(pm.Categorical):
    """
    Hidden Markov Model States
    Parameters
    ----------
    P1 : tensor
        probability to remain in state 1
    P2 : tensor
        probability to move from state 2 to state 1

    """

    def __init__(self, PA=None, P1=None, P2=None,
                 *args, **kwargs):
        super(pm.Categorical, self).__init__(*args, **kwargs)
        self.PA = PA
        self.P1 = P1
        self.P2 = P2
        self.k = 2 # two state model
        self.mean = 0.
        self.mode = tt.cast(0,dtype='int64')

    def logp(self, x):
        PA = self.PA
        P1 = self.P1
        P2 = self.P2

        # now we need to create an array with probabilities
        # so that for x=A: PA=P1, PB=(1-P1)
        # and for x=B: PA=P2, PB=(1-P2)
        PT = tt.stack((P1,P2))

        P = PT[x[:-1]]

        x_i = x[1:]

        ou_like = pm.Categorical.dist(P, shape=(N_chain-1,2)).logp(x_i)
        return pm.Categorical.dist(PA).logp(x[0]) + tt.sum(ou_like)

我的 HMMStates 不能真正调用 pm.Categorical 超级初始化，所以我调用 pm.Categorical 的超级 class 即 pm.Discrete。这个技巧使它通过了 BinaryGibbsMetropolis 和 CategoricalGibbsMetropolis 的测试。

如果您有兴趣实现 2 状态和多状态 HMM，我将实现所有这些案例 here。

我可以使用什么 pymc3 Monte-Carlo 步进器进行自定义分类分布？

What pymc3 Monte-Carlo stepper can I use for a custom categorical distribution?

bayesian

hidden-markov-models

theano

pymc3