Как сделать марковскую модель с дискретным состоянием с помощью pymc?

Я пытаюсь понять, как правильно сделать модель цепи Маркова с дискретным состоянием с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])

Выборка из модели показывает, что маргинальные состояния состояний должны быть такими (по сравнению с моделью, созданной с помощью пакета BNT Кевина Мерфи в Matlab)

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 для более общих динамических байесовских сетей.