Какой шагомер pymc3 Monte-Carlo можно использовать для пользовательского категориального распространения?

Question

Я работаю над реализацией скрытых Марков-цепей в pymc3. Я довольно успешно реализовал скрытые состояния. Ниже я показываю простой 2-состояния марковской цепью:

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) 

Я теперь определив класс, который описывает состояния. В качестве ввода мне нужно знать, что вероятность P1 переместится из состояния 0 в состояние 1, а P2 - от 1-> 0. Мне также нужно знать вероятность PA для первого состояния будучи 0.

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, чтобы избежать ошибки testvalue. Вот как использовать класс в модели, которая проверяет, работает ли она. В этом случае состояние не скрыто, но наблюдается.

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) 

выход воспроизводит данные красиво. В следующей модели я покажу проблему. Здесь я непосредственно не наблюдаю состояния, но состояние с некоторым гауссовским шумом добавлено (таким образом, скрытое состояние). Если вы запустите модель с помощью Stepper Metropolis, она обрушится с ошибкой индекса, и я вернусь к проблемам, связанным с использованием Stepper Metropolis в категориальных дистрибутивах. К сожалению, единственным шагером, который применим к моему классу, является степлер CategoricalGibbsMetropolis, но он отказывается работать с моим классом, поскольку он явно не является категориальным распределением.

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, либо все 1s.

Как я могу сказать ElemwiseCategorial, чтобы назначить вектор состояний из 1s и 0s, или, как альтернатива, как получить CategorialGibbsMetropolis, чтобы признать мое распределение как категориальное. Это должна быть общая проблема с пользовательскими дистрибутивами.

Answer 1

Поскольку я ни от кого не слышал по моему вопросу, я сам ответил сам. Трюк, который я использовал, был предложен Крисом Фоннебеком на pymc3 github, где я открыл проблему. Он предложил подкласс pm. Категориальный.

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, который pm.Discrete. Этот трюк заставляет пройти тест BinaryGibbsMetropolis и CategoricalGibbsMetropolis.

Если вы заинтересованы в внедрении HMM с состоянием состояния из нескольких состояний и нескольких состояний, я буду выполнять все эти случаи here.