如何用pymc构建离散状态马尔可夫模型？

Question

我正在尝试弄清楚如何用 pymc 正确地创建一个离散状态的马尔可夫链模型。

举个例子（可以在 nbviewer 查看），我们来创建一个长度为 T=10 的链，其中马尔可夫状态是二进制的，初始状态分布是 [0.2, 0.8]，在状态 1 时切换状态的概率是 0.01，而在状态 2 时是 0.5。

import numpy as np
import pymc as pm
T = 10
prior0 = [0.2, 0.8]
transMat = [[0.99, 0.01], [0.5, 0.5]]

为了构建这个模型，我创建了一个状态变量的数组和一个依赖于状态变量的转移概率数组（使用 pymc.Index 函数）。

states = np.empty(T, dtype=object)
states[0] = pm.Categorical('state_0', prior0)
transPs = np.empty(T, dtype=object)
transPs[0] = pm.Index('trans_0', transMat, states[0])

for i in range(1, T):
    states[i] = pm.Categorical('state_%i' % i, transPs[i-1])
    transPs[i] = pm.Index('trans_%i' %i, transMat, states[i])

对模型进行采样显示，状态的边际分布是正确的（与在 Matlab 中使用 Kevin Murphy 的 BNT 包构建的模型相比）。

model = pm.MCMC([states, transPs])
model.sample(10000, 5000)
[np.mean(model.trace('state_%i' %i)[:]) for i in range(T)]    

输出结果是：

[-----------------100%-----------------] 10000 of 10000 complete in 7.5 sec

[0.80020000000000002,
 0.39839999999999998,
 0.20319999999999999,
 0.1118,
 0.064199999999999993,
 0.044600000000000001,
 0.033000000000000002,
 0.026200000000000001,
 0.024199999999999999,
 0.023800000000000002]

我的问题是——这似乎不是用 pymc 构建马尔可夫链的最优雅的方法。有没有更简洁的方法，不需要使用确定性函数的数组？

我的目标是编写一个基于 pymc 的包，用于更通用的动态贝叶斯网络。

Answer 1

如果你想了解你的马尔可夫链在长时间运行后的表现，discreteMarkovChain 这个包可能会对你有帮助。里面的例子展示了如何通过定义一个转移函数来构建一个离散状态的马尔可夫链。这个转移函数可以告诉你在每个状态下，下一步可以到达哪些状态以及它们的概率。你也可以用这个函数来模拟这个过程。

Answer 2

据我所知，你需要把每个时间点的分布编码成前一个时间点的确定性函数，因为这就是它的本质——参数中没有随机性，因为你在问题设置中已经定义好了它们。不过，我觉得你的问题可能是想找一种更直观的方式来表示这个模型。 一种替代的方法是直接把时间点之间的转换编码成前一个时间点的函数。

from pymc import Bernoulli, MCMC

def generate_timesteps(N,p_init,p_trans):
    timesteps=np.empty(N,dtype=object)
    # A success denotes being in state 2, a failure being in state 1
    timesteps[0]=Bernoulli('T0',p_init)
    for i in xrange(1,N):
        # probability of being in state 1 at time step `i` given time step `i-1`
        p_i = p_trans[1]*timesteps[i-1]+p_trans[0]*(1-timesteps[i-1])
        timesteps[i] = Bernoulli('T%d'%i,p_i)
    return timesteps

timesteps = generate_timesteps(10,0.8,[0.001,0.5])
model = MCMC(timesteps)
model.sample(10000) # no burn in necessary since we're sampling directly from the distribution
[np.mean( model.trace(t).gettrace() ) for t in timesteps]