Можно ли использовать операторы if в детерминированной функции PyMC?

Question

После прочтения вероятностного программирования и байесовских методов для хакеров Кэма Дэвидсона-Пилона я решил попробовать свои силы в задаче изучения скрытой марковской модели (HMM) с PyMC. Пока что код не взаимодействует, но из-за устранения неполадок я чувствую, что сузил источник проблемы.

Разбивая код на более мелкие порции и сосредотачиваясь на начальной вероятности и вероятностях выбросов при t=0, я могу узнать параметры излучения / наблюдения одного состояния в момент времени t=0. Однако, как только я добавляю в другое состояние (всего два состояния), результаты обучения параметров идентичны (и неверны) независимо от ввода данных. Итак, я чувствую, что я сделал что-то не так в @pm.deterministic часть кода, которая не позволяет мне пробовать из Init начальная функция вероятности.

С помощью этой части кода я стремлюсь узнать начальную вероятность p_bern и вероятности выбросов p_0 а также p_1 соответствует состояниям 0 и 1 соответственно. Эмиссия зависит от состояния, которое я пытаюсь выразить своим @pm.deterministic функция. Могу ли я использовать оператор if в этой детерминированной функции? Кажется, это корень проблемы.

# This code is to test the ability to discern between two states with emissions

import numpy as np
import pymc as pm
from matplotlib import pyplot as plt

N = 1000
state = np.zeros(N)
data = np.zeros(shape=N)

# Generate data
for i in range(N):
    state[i] = pm.rbernoulli(p=0.3)
for i in range(N):
    if state[i]==0:
        data[i] = pm.rbernoulli(p=0.4)
    elif state[i]==1:
        data[i] = pm.rbernoulli(p=0.8)

# Prior on probabilities
p_bern = pm.Uniform("p_S", 0., 1.)
p_0 = pm.Uniform("p_0", 0., 1.)
p_1 = pm.Uniform("p_1", 0., 1.)

Init = pm.Bernoulli("Init", p=p_bern) # Bernoulli node

@pm.deterministic
def p_T(Init=Init, p_0=p_0, p_1=p_1, p_bern=p_bern):
    if Init==0:
        return p_0
    elif Init==1:
        return p_1

obs = pm.Bernoulli("obs", p=p_T, value=data, observed=True)
model = pm.Model([obs, p_bern, p_0, p_1])
mcmc = pm.MCMC(model)
mcmc.sample(20000, 10000)
pm.Matplot.plot(mcmc)

Я уже пытался следующее безрезультатно:

1. Использовать @pm.potential декоратор для создания совместного распространения
2. Изменение места расположения моего Init местоположение (вы можете увидеть мой комментарий в коде, где я не уверен, где его разместить)
3. Использовать @pm.stochastic похоже на это

Редактировать: Согласно предложению Криса, я переместил узел Бернулли за пределы детерминизма. Я также обновил код для упрощения модели (наблюдение Бернулли вместо многочлена) для упрощения поиска неисправностей.

Спасибо за ваше время и внимание. Любые отзывы тепло принимаются. Кроме того, если мне не хватает какой-либо информации, пожалуйста, дайте мне знать!

Answer 1

Я бы убрал эту стохастичность из детерминизма. Предполагается, что значение детерминированного узла полностью определяется значением его родителей. Наличие случайной переменной, скрытой в узле, нарушает это.

Зачем создавать точку Бернулли и передавать ее в качестве аргумента детерминистическому?

Answer 2

Основываясь на предоставленной вами обновленной информации, вот код, который работает:

import numpy as np
import pymc as pm
from matplotlib import pyplot as plt

N = 1000
state = np.zeros(N)
data = np.zeros(shape=N)

# Generate data
state = pm.rbernoulli(p=0.3, size=N)
data = [int(pm.rbernoulli(0.8*s or 0.4)) for s in state]

# Prior on probabilities
p_S = pm.Uniform("p_S", 0., 1.)
p_0 = pm.Uniform("p_0", 0., 1.)
p_1 = pm.Uniform("p_1", 0., 1.)

# Use values of Init as indices to probabilities
Init = pm.Bernoulli("Init", p=p_S, size=N) # Bernoulli node

p_T = pm.Lambda('p_T', lambda p_0=p_0, p_1=p_1, i=Init: np.array([p_0, p_1])[i.astype(int)])

obs = pm.Bernoulli("obs", p=p_T, value=data, observed=True)
model = pm.MCMC(locals())
model.sample(20000, 10000)
model.summary()

Обратите внимание, что на этапе генерации данных я использую состояния для индексации соответствующей истинной вероятности. По сути, я делаю то же самое в спецификации p_T, Кажется, он работает достаточно хорошо, но учтите, что в зависимости от того, где что-то инициализируется, два значения p_0 а также p_1 может в конечном итоге соответствовать одному из истинных значений (ничто не может ограничить одно быть больше другого). Итак, ценность p_S может закончиться как дополнение к вероятности истинного состояния.