Почему Pymc3 ADVI хуже, чем MCMC в этом примере логистической регрессии? - PullRequest
0 голосов
/ 28 сентября 2018

Мне известны математические различия между ADVI / MCMC, но я пытаюсь понять практические последствия использования того или другого.Я запускаю очень простой пример логистической регрессии для данных, которые я создал таким образом:

import pandas as pd
import pymc3 as pm
import matplotlib.pyplot as plt
import numpy as np

def logistic(x, b, noise=None):
    L = x.T.dot(b)
    if noise is not None:
        L = L+noise
    return 1/(1+np.exp(-L))

x1 = np.linspace(-10., 10, 10000)
x2 = np.linspace(0., 20, 10000)
bias = np.ones(len(x1))
X = np.vstack([x1,x2,bias]) # Add intercept
B =  [-10., 2., 1.] # Sigmoid params for X + intercept

# Noisy mean
pnoisy = logistic(X, B, noise=np.random.normal(loc=0., scale=0., size=len(x1)))
# dichotomize pnoisy -- sample 0/1 with probability pnoisy
y = np.random.binomial(1., pnoisy)

И я запускаю ADVI следующим образом:

with pm.Model() as model: 
    # Define priors
    intercept = pm.Normal('Intercept', 0, sd=10)
    x1_coef = pm.Normal('x1', 0, sd=10)
    x2_coef = pm.Normal('x2', 0, sd=10)

    # Define likelihood
    likelihood = pm.Bernoulli('y',                  
           pm.math.sigmoid(intercept+x1_coef*X[0]+x2_coef*X[1]),
                          observed=y)
    approx = pm.fit(90000, method='advi')

К сожалению, сколько бы я ниПри увеличении выборки ADVI не может восстановить исходные бета-версии, которые я определил [-10., 2., 1.], а MCMC работает нормально (как показано ниже)

enter image description here

Спасибо за помощь!

1 Ответ

0 голосов
/ 29 сентября 2018

Это интересный вопрос!Значение по умолчанию 'advi' в PyMC3 - это средний вариационный вывод поля, который не очень хорошо справляется с захватом корреляций.Оказывается, что настроенная вами модель имеет интересную структуру корреляции, которую можно увидеть с помощью этого:

import arviz as az

az.plot_pair(trace, figsize=(5, 5))

correlated samples

PyMC3 имеет встроенную- проверка сходимости - выполнение оптимизации на длинные или слишком короткие может привести к забавным результатам:

from pymc3.variational.callbacks import CheckParametersConvergence

with model:
    fit = pm.fit(100_000, method='advi', callbacks=[CheckParametersConvergence()])

draws = fit.sample(2_000)

Это останавливается после примерно 60 000 итераций для меня.Теперь мы можем проверить корреляции и увидеть, что, как и ожидалось, ADVI соответствует гауссианам, выровненным по оси:

az.plot_pair(draws, figsize=(5, 5))

another correlation image

Наконец, мы можем сравнитьподходит от NUTS и (среднее поле) ADVI:

az.plot_forest([draws, trace])

forest plot

Обратите внимание, что ADVI недооценивает дисперсию, но достаточно близко для среднего значения каждого параметра,Кроме того, вы можете установить method='fullrank_advi', чтобы фиксировать корреляции, которые вы видите немного лучше.

(примечание: arviz скоро станет библиотекой черчения для PyMC3)

...