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

Выборка из модели показывает, что маргинальные состояния состояний должны быть такими, какими они должны быть (по сравнению с моделью, построенной с помощью пакета 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 для более общих динамических байесовских сетей.

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]