Изучение дискретных параметров HMM в PyMC

Question

Я пытаюсь узнать параметры простого дискретного HMM, используя PyMC. Я моделирую дождливо-солнечную модель со страницы Wiki на HMM. Модель выглядит следующим образом:

Я использую следующие приоры.

theta_start_state ~ beta(20,10)
theta_transition_rainy ~beta(8,2)
theta_transition_sunny ~beta(2,8)
theta_emission_rainy ~ Dirichlet(3,4,3)
theta_emission_sunny ~ Dirichlet(10,6,4)

Первоначально я использую эту настройку для создания обучающего набора следующим образом.

## Some not so informative priors!
# Prior on start state
theta_start_state = pm.Beta('theta_start_state',12,8)

# Prior on transition from rainy
theta_transition_rainy = pm.Beta('transition_rainy',8,2)

# Prior on transition from sunny
theta_transition_sunny = pm.Beta('transition_sunny',2,8)

# Prior on emission from rainy
theta_emission_rainy = pm.Dirichlet('emission_rainy',[3,4,3])

# Prior on emission from sunny
theta_emission_sunny = pm.Dirichlet('emission_sunny',[10,6,4])

# Start state
x_train_0 = pm.Categorical('x_0',[theta_start_state, 1-theta_start_state])

N = 100

# Create a train set for hidden states
x_train = np.empty(N, dtype=object)

# Creating a train set of observations
y_train = np.empty(N, dtype=object)

x_train[0] = x_train_0

for i in xrange(1, N):
    if x_train[i-1].value==0:
        x_train[i] = pm.Categorical('x_train_%d'%i,[theta_transition_rainy, 1- theta_transition_rainy])
    else:
        x_train[i] = pm.Categorical('x_train_%d'%i,[theta_transition_sunny, 1- theta_transition_sunny])


for i in xrange(0,N):
    if x_train[i].value == 0:
        # Rain
        y_train[i] = pm.Categorical('y_train_%d' %i, theta_emission_rainy)
    else:
        y_train[i] = pm.Categorical('y_train_%d' %i, theta_emission_sunny)

Тем не менее, я не могу понять, как узнать эти параметры с помощью PyMC. Я сделал старт следующим образом.

@pm.observed
def y(x=x_train, value =y_train):    
    N = len(x)    
    out = np.empty(N, dtype=object)
    for i in xrange(0,N):
        if x[i].value == 0:
            # Rain
            out[i] = pm.Categorical('y_%d' %i, theta_emission_rainy)
        else:
            out[i] = pm.Categorical('y_%d' %i, theta_emission_sunny)
    return out

Полный блокнот, содержащий этот код, можно найти здесь.

За исключением: суть, содержащая HMM-код для гауссиана, действительно трудно понять! (не задокументировано)

Обновить

Основываясь на ответах ниже, я попытался изменить свой код следующим образом:

@pm.stochastic(observed=True)
def y(value=y_train, hidden_states = x_train):
    def logp(value, hidden_states):
        logprob = 0 
        for i in xrange(0,len(hidden_states)):
            if hidden_states[i].value == 0:
                # Rain
                logprob = logprob + pm.categorical_like(value[i], theta_emission_rainy)
            else:
                # Sunny
                logprob = logprob + pm.categorical_like(value[i], theta_emission_sunny)
        return logprob 

Следующим шагом будет создание модели и запуск алгоритма MCMC. Однако отредактированный выше код также не будет работать. Это дает ZeroProbability error,

Я не уверен, правильно ли я истолковал ответы.

Answer 1

Просто некоторые мысли по этому поводу:

• Санни и Рейни - взаимоисключающие и исчерпывающие скрытые состояния. Почему бы вам не кодировать их как одну категориальную переменную погоды, которая может принимать одно из двух значений (кодирование для солнечного, дождливого)?
• В вашей функции правдоподобия вы, кажется, наблюдаете Дождь / Солнечный. То, как я вижу это на графике вашей модели, похоже, что это скрытые, а не наблюдаемые переменные (это были бы "прогулка", "магазин" и "чистота")

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

Учусь

Если вы хотите узнать параметры для модели, вы можете рассмотреть возможность перехода на PyMC3, который лучше подходит для автоматического вычисления градиентов для вашей функции logp. Но в этом случае (поскольку вы выбрали сопряженные априорные значения) это не является действительно необходимым. Если вы не знаете, что такое Conjugate Priors, или вам нужен обзор, попросите Wikipedia для получения списка Conjugate Priors, у вас есть отличная статья на эту тему.

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

Если вас не интересует маргинальное апостериорное распределение, а, скорее, поиск совместных параметров MAP, вы можете рассмотреть возможность использования обучения с максимизацией ожидания (EM) или имитацией отжига. Оба должны работать достаточно хорошо в рамках MCMC Framework.

Для EM Learning просто повторите эти шаги до сходимости:

• E (Ожидание) Шаг: Запустите на некоторое время цепочку MCMC и соберите образцы
• Шаг M (максимизация): Обновите гиперпараметры для ваших бета-версий и дирихле-приоров, как если бы эти образцы были фактическими наблюдениями. Посмотрите дистрибутивы Beta и Dirichlet, если вы не знаете, как это сделать.

Я бы использовал небольшой коэффициент скорости обучения, чтобы вы не попали в первый локальный оптимум (теперь мы приближаемся к моделируемому отжигу): вместо того, чтобы рассматривать N выборок, которые вы генерировали из цепочки MCMC, как фактические наблюдения, трактуйте их как K наблюдения (для значения K << N) путем масштабирования обновлений гиперпараметров вниз на коэффициент скорости обучения K/N.

Answer 2

Первое, что бросается в глаза - это возвращаемое значение вашей вероятности. PyMC ожидает скалярное возвращаемое значение, а не список / массив. Вам необходимо сложить массив перед его возвратом.

Кроме того, когда вы используете Dirichlet в качестве априора для Категориального, PyMC обнаруживает это и заполняет последнюю вероятность. Вот как я бы закодировал ваш x_train/y_train петли:

p = []
for i in xrange(1, N):
    # This will return the first variable if prev=0, and the second otherwise
    p.append(pm.Lambda('p_%i' % i, lambda prev=x_train[i-1]: (theta_transition_rainy, theta_transition_sunny)[bool(prev)]))
    x_train[i] = pm.Categorical('x_train_%i' % i, p[-1])

Итак, вы берете соответствующие вероятности с помощью лямбды и используете ее в качестве аргумента для категориального.