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


python markov-chains pymc
person shpigi    schedule 25.03.2014    source источник

Ответы (2)


arrow_upward
2
arrow_downward

Насколько я знаю, вы должны кодировать распределение каждого временного шага как детерминированную функцию предыдущего временного шага, потому что это то, что есть — в параметрах нет случайности, потому что вы определили их в постановке задачи. Однако я думаю, что ваш вопрос, возможно, был больше связан с поиском более интуитивно понятного способа представления модели. Одним из альтернативных способов было бы прямое кодирование переходов временного шага в зависимости от предыдущего временного шага.

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]
person Austen    schedule 25.03.2015

arrow_upward
2
arrow_downward

Если вы хотите посмотреть на долгосрочное поведение вашей цепи Маркова, пакет discreteMarkovChain может быть полезный. В примерах показаны некоторые идеи построения цепи Маркова с дискретными состояниями путем определения функции перехода, которая сообщает вам для каждого состояния достижимые состояния на следующем шаге и их вероятности. Вы можете использовать ту же функцию для имитации процесса.

person Forzaa    schedule 09.04.2016