Код PyMC дает необычные результаты

Question

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

Модель выглядит следующим образом:

Y_i = Bernoulli(p_i) 
logit(p_i) = b0 + b1*x1 + b2*x2 

Реализация заключается в следующем:

import pymc as pm 
import numpy as np 

b0 = pm.Normal('b0', 0., 1e-6, value=0.) 
b1 = pm.Normal('b1', 0., 1e-6, value=0.) 
b2 = pm.Normal('b2', 0., 1e-6, value=0.) 

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32]) 
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15]) 
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0]) 

@pm.deterministic 
def p(b0=b0, b1=b1, b2=b2, x1=x1, x2=x2): 
    return np.exp(b0 + b1*x1 + b2*x2)/(1.0 + np.exp(b0 + b1*x1 + b2*x2)) 

obs = pm.Bernoulli('obs', p=p, value=value, observed=True) 

m = pm.MCMC([obs, b0, b1, b2]) 

Когда я сэмплирую m.sample(500000, 200000, 50) и построить полученное распределение задний, я получаю это:

Во второй попытке получить лучшие результаты, я использовал @ pm.observed:

import pymc as pm 
import numpy as np 

b0 = pm.Normal('b0', 0., 1e-6, value=0.) 
b1 = pm.Normal('b1', 0., 1e-6, value=0.) 
b2 = pm.Normal('b2', 0., 1e-6, value=0.) 

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32]) 
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15]) 
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0]) 

p = b0 + b1*x1 + b2*x2 

@pm.observed 
def obs(p=p, value=value): 
    return pm.bernoulli_like(value, pm.invlogit(p)) 

m = pm.MCMC([obs, p, b0, b1, b2]) 

Но она также производит высокие автокорреляции.

Я увеличил размер выборки без особого успеха. Что мне не хватает?

Answer 1

У вас определенно есть проблемы с конвергенцией. Частично проблема заключается в том, что у вас очень маленький размер выборки (n = 9), но вы пытаетесь установить 3 параметра. Я получаю немного лучше пробег по блочно-обновлению параметров, но перехватывать все еще плохо работает:

beta = pm.Normal('beta', 0., 0.001, value=[0.]*3) 

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32]) 
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15]) 
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0]) 

p = pm.Lambda('p', lambda b=beta: pm.invlogit(b[0] + b[1]*x1 + b[2]*x2)) 

obs = pm.Bernoulli('obs', p=p, value=value, observed=True) 

(Кстати, когда вы логит-преобразования ваших параметров, вам не нужны супер-диффузные априорные на вашем бета - меньше 0,001 - довольно бесплодный). Вещи улучшаются дальше, используя адаптивный Метрополис пробник на бетах:

M.use_step_method(AdaptiveMetropolis, M.beta) 

Это приводит к конвергенции, но, как вы можете видеть, что нет никакой информации для информирования перед: