Байесовский лассо с использованием PyMC3

Я пытаюсь воспроизвести результаты этого tutorial (см. Регрессию LASSO) на PyMC3. Как прокомментировал этот reddit thread, смешивание для первых двух коэффициентов не было хорошим, потому что переменные коррелированы.Байесовский лассо с использованием PyMC3

Я попытался реализовать его в PyMC3, но он не работал должным образом при использовании гамильтоновых пробоотборников. Я мог только заставить его работать с пробником Metropolis, который достигает того же результата, что и PyMC2.

Я не знаю, связано ли это с тем, что лапласиан достиг максимума (разрывная производная при 0), но он отлично работал с гауссовыми приурами. Я пытался с инициализацией MAP или без нее, и результат всегда один и тот же.

Вот мой код:

from pymc import * 
from scipy.stats import norm 
import pylab as plt 

# Same model as the tutorial 
n = 1000 

x1 = norm.rvs(0, 1, size=n) 
x2 = -x1 + norm.rvs(0, 10**-3, size=n) 
x3 = norm.rvs(0, 1, size=n) 

y = 10 * x1 + 10 * x2 + 0.1 * x3 

with Model() as model: 
    # Laplacian prior only works with Metropolis sampler 
    coef1 = Laplace('x1', 0, b=1/sqrt(2)) 
    coef2 = Laplace('x2', 0, b=1/sqrt(2)) 
    coef3 = Laplace('x3', 0, b=1/sqrt(2)) 

    # Gaussian prior works with NUTS sampler 
    #coef1 = Normal('x1', mu = 0, sd = 1) 
    #coef2 = Normal('x2', mu = 0, sd = 1) 
    #coef3 = Normal('x3', mu = 0, sd = 1) 

    likelihood = Normal('y', mu= coef1 * x1 + coef2 * x2 + coef3 * x3, tau = 1, observed=y) 

    #step = Metropolis() # Works just like PyMC2 
    start = find_MAP() # Doesn't help 
    step = NUTS(state = start) # Doesn't work 
    trace = sample(10000, step, start = start, progressbar=True) 

plt.figure(figsize=(7, 7)) 
traceplot(trace) 
plt.tight_layout() 

autocorrplot(trace) 
summary(trace)

Здесь ошибка я получаю:

PositiveDefiniteError: Simple check failed. Diagonal contains negatives

Я делаю что-то неправильно или орехи пробник не должен работать на подобных случаях?

источник

2015-02-03 Pedro Tabacof

Я пошел вперед и переписал модель в Стэн. Он не сломается, но для запуска NUTS требуется очень много времени. Я подозреваю, что, поскольку модель не учитывает тот факт, что x1 ~ = -x2, есть просто невероятное количество отказов. Поэтому, чтобы ответить на ваш вопрос: я не считаю, что NUTS должны работать над этим делом. –

«Лапласиан» делает это Лассо? –

@ O.rka Lasso - это решение максимальной линейки с линейной моделью с приоритетом и лапласией. Если вместо MAP вы делаете байесовский вывод, это не совсем лассо, но это связано. См. Эту статью для получения дополнительной информации: http://www.stat.ufl.edu/archived/casella/Papers/Lasso.pdf –

Whyking из reddit thread дал предложение использовать MAP как масштабирование вместо состояния, и это действительно сработало чудеса.

Адрес notebook с результатами и обновленными code.

источник

2015-02-03 22:25:28

Байесовский лассо с использованием PyMC3

ответ

Смежные вопросы