PyMC2 и PyMC3 дают разные результаты ...?

Question

Я пытаюсь получить простую модель PyMC2, работающую в PyMC3. Я получил модель для запуска, но модели дают очень разные оценки MAP для переменных. Вот моя модель PyMC2:

import pymc 
theta = pymc.Normal('theta', 0, .88) 

X1 = pymc.Bernoulli('X2', p=pymc.Lambda('a', lambda theta=theta:1./(1+np.exp(-(theta-(-0.75))))), value=[1],observed=True) 
X2 = pymc.Bernoulli('X3', p=pymc.Lambda('b', lambda theta=theta:1./(1+np.exp(-(theta-0)))), value=[1],observed=True) 

model = pymc.Model([theta, X1, X2]) 
mcmc = pymc.MCMC(model) 
mcmc.sample(iter=25000, burn=5000) 
trace = (mcmc.trace('theta')[:]) 
print "\nThe MAP value for theta is", trace.sum()/len(trace) 

Это похоже на работу, как ожидалось. У меня были всевозможные проблемы с выяснением того, как использовать эквивалент объекта pymc.Lambda в PyMC3. В конце концов я столкнулся с детерминированным объектом. Ниже мой код:

import pymc3 

with pymc3.Model() as model: 

    theta = pymc3.Normal('theta', 0, 0.88) 
    X1 = pymc3.Bernoulli('X1', p=pymc3.Deterministic('b', 1./(1+np.exp(-(theta-(-0.75))))), observed=[1]) 
    X2 = pymc3.Bernoulli('X2', p=pymc3.Deterministic('c', 1./(1+np.exp(-(theta-(0))))), observed=[1]) 

    start=pymc3.find_MAP() 
    step=pymc3.NUTS(state=start) 
    trace = pymc3.sample(20000, step, njobs=1, progressbar=True) 

pymc3.traceplot(trace) 

Проблема у меня в том, что моя оценка MAP для тета с использованием PyMC2 составляет ~ 0,68 (правильно), в то время как оценка PyMC3 дает составляет ~ 0,26 (неправильно). Я подозреваю, что это имеет какое-то отношение к тому, как я определяю детерминированную функцию. PyMC3 не позволит мне использовать лямбда-функцию, поэтому мне просто нужно написать выражение в строке. Когда я пытаюсь использовать лямбда-тета = тета: ... Я получаю эту ошибку:

AsTensorError: ('Cannot convert <function <lambda> at 0x157323e60> to TensorType', <type 'function'>) 

что-то делать с Теано ?? Любые предложения будут ценны!

Answer 1

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

• find_MAP дал неправильный ответ
• орехов пробоотборник дал неверному ответу
• Метрополис пробнику дал правильный ответ, яй!

Я читал где-то еще, что сэмплер NUTS не работает с детерминированным. Я не знаю почему. Может быть, так обстоит дело с find_MAP? Но пока я буду придерживаться Метрополиса.

Answer 2

Это работает, когда вы используете anano tensor вместо функции numpy в вашем детерминированном режиме.

import pymc3 
import theano.tensor as tt 

with pymc3.Model() as model: 

    theta = pymc3.Normal('theta', 0, 0.88) 
    X1 = pymc3.Bernoulli('X1', p=pymc3.Deterministic('b', 1./(1+tt.exp(-(theta-(-0.75))))), observed=[1]) 
    X2 = pymc3.Bernoulli('X2', p=pymc3.Deterministic('c', 1./(1+tt.exp(-(theta-(0))))), observed=[1]) 

    start=pymc3.find_MAP() 
    step=pymc3.NUTS(state=start) 
    trace = pymc3.sample(20000, step, njobs=1, progressbar=True) 

print "\nThe MAP value for theta is", np.median(trace['theta']) 

pymc3.traceplot(trace); 

Вот результат:

Answer 3

Кроме того, NUTS не обрабатывает дискретные переменные. Если вы хотите использовать NUTS, вы должны разделить пробоотборников:

step1 = pymc3.NUTS([theta]) 
step2 = pymc3.BinaryMetropolis([X1,X2]) 

trace = pymc3.sample(10000, [step1, step2], start) 

EDIT: Пропущенный, что «б» и «с» были определены инлайн. Удалены из функции функции NUTS

Answer 4

Значение MAP не определяется как среднее для распределения, а как его максимальное значение. С pymc2 вы можете найти его с:

M = pymc.MAP(model) 
M.fit() 
theta.value 

который возвращает array(0.6253614422469552)

Это согласуется с МАП, которые вы найдете с find_MAP в pymc3, которую вы называете start:

{'theta': array(0.6253614811102668)} 

Вопрос о который является лучшим сэмплером, отличается от него и не зависит от вычисления MAP.Расчет MAP является оптимизацией. См.: https://pymc-devs.github.io/pymc/modelfitting.html#maximum-a-posteriori-estimates для pymc2.