PyMC3 Байесовское линейное регрессионное прогнозирование с использованием sklearn.datasets

Я пытаюсь реализовать модели Байесовской линейной регрессии, используя PyMC3 с REAL DATA (т.е. не из линейной функции + гауссовский шум) из наборов данных в sklearn.datasets. Я выбрал набор регрессионных данных с наименьшим количеством атрибутов (т.е. load_diabetes()), форма которого (442, 10); то есть 442 samples и 10 attributes.

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

Я пробовал: Создание прогнозов из предполагаемых параметров в pymc3 а также http://pymc-devs.github.io/pymc3/posterior_predictive/, но моя модель либо чрезвычайно ужасна при прогнозировании, либо я делаю это неправильно.

Если я действительно правильно делаю прогноз (чего я, вероятно, нет), то кто-нибудь может помочь мне оптимизировать мою модель. Я не знаю, если меньше mean squared error, absolute error, или что-то в этом роде работает в байесовских рамках. В идеале я хотел бы получить массив number_of_rows = количество строк в наборе тестов атрибутов/данных X_te и количество столбцов, которые будут выборками из заднего распределения.

import pymc3 as pm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from scipy import stats, optimize
from sklearn.datasets import load_diabetes
from sklearn.cross_validation import train_test_split
from theano import shared

np.random.seed(9)

%matplotlib inline

#Load the Data
diabetes_data = load_diabetes()
X, y_ = diabetes_data.data, diabetes_data.target

#Split Data
X_tr, X_te, y_tr, y_te = train_test_split(X,y_,test_size=0.25, random_state=0)

#Shapes
X.shape, y_.shape, X_tr.shape, X_te.shape
#((442, 10), (442,), (331, 10), (111, 10))

#Preprocess data for Modeling
shA_X = shared(X_tr)

#Generate Model
linear_model = pm.Model()

with linear_model: 
    # Priors for unknown model parameters    
    alpha = pm.Normal("alpha", mu=0,sd=10)
    betas = pm.Normal("betas", mu=0,#X_tr.mean(), 
                               sd=10, 
                               shape=X.shape[1])
    sigma = pm.HalfNormal("sigma", sd=1)

    # Expected value of outcome
    mu = alpha + np.array([betas[j]*shA_X[:,j] for j in range(X.shape[1])]).sum()

    # Likelihood (sampling distribution of observations)
    likelihood = pm.Normal("likelihood", mu=mu, sd=sigma, observed=y_tr)

    # Obtain starting values via Maximum A Posteriori Estimate
    map_estimate = pm.find_MAP(model=linear_model, fmin=optimize.fmin_powell)

    # Instantiate Sampler
    step = pm.NUTS(scaling=map_estimate)

    # MCMC
    trace = pm.sample(1000, step, start=map_estimate, progressbar=True, njobs=1)


#Traceplot
pm.traceplot(trace)

введите описание изображения здесь

# Prediction
shA_X.set_value(X_te)
ppc = pm.sample_ppc(trace, model=linear_model, samples=1000)

#What the shape of this? 
list(ppc.items())[0][1].shape #(1000, 111) it looks like 1000 posterior samples for the 111 test samples (X_te) I gave it

#Looks like I need to transpose it to get `X_te` samples on rows and posterior distribution samples on cols

for idx in [0,1,2,3,4,5]:
    predicted_yi = list(ppc.items())[0][1].T[idx].mean()
    actual_yi = y_te[idx]
    print(predicted_yi, actual_yi)
# 158.646772735 321.0
# 160.054730647 215.0
# 149.457889418 127.0
# 139.875149489 64.0
# 146.75090354 175.0
# 156.124314452 275.0

Ответы

Ответ 1

Я думаю, что одна из проблем с вашей моделью заключается в том, что ваши данные имеют очень разные масштабы, у вас есть диапазон ~ 0.3 для ваших "Xs" и ~ 300 для вашего "Ys". Следовательно, вы должны ожидать больших склонов (и сигмы), которые ваши предвидения задают. Один логический вариант - настроить ваши приоритеты, как в следующем примере.

#Generate Model
linear_model = pm.Model()

with linear_model: 
    # Priors for unknown model parameters    
    alpha = pm.Normal("alpha", mu=y_tr.mean(),sd=10)
    betas = pm.Normal("betas", mu=0, sd=1000, shape=X.shape[1])
    sigma = pm.HalfNormal("sigma", sd=100) # you could also try with a HalfCauchy that has longer/fatter tails
    mu = alpha + pm.dot(betas, X_tr.T)
    likelihood = pm.Normal("likelihood", mu=mu, sd=sigma, observed=y_tr)
    step = pm.NUTS()
    trace = pm.sample(1000, step)

chain = trace[100:]
pm.traceplot(chain);

введите описание изображения здесь

Задние прогностические проверки показывают, что у вас более или менее разумная модель.

sns.kdeplot(y_tr, alpha=0.5, lw=4, c='b')
for i in range(100):
    sns.kdeplot(ppc['likelihood'][i], alpha=0.1, c='g')

введите описание изображения здесь

Другой вариант заключается в том, чтобы поместить данные в один и тот же масштаб, стандартизируя его, сделав это, вы получите, что наклон должен быть около + -1, и в общем случае вы можете использовать тот же самый диффуз для любых данных (что-то полезное, если только у вас есть информативные приоритеты, которые вы можете использовать). Фактически, многие люди рекомендуют эту практику для обобщенных линейных моделей. Вы можете больше узнать об этом в книге анализ балийских данных или Статистическое переосмысление

Если вы хотите предсказать значения, у вас есть несколько вариантов, нужно использовать среднее значение выводимых параметров, например:

alpha_pred = chain['alpha'].mean()
betas_pred = chain['betas'].mean(axis=0)

y_pred = alpha_pred + np.dot(betas_pred, X_tr.T)

Другой вариант - использовать pm.sample_ppc для получения образцов прогнозируемых значений, которые учитывают неопределенность в ваших оценках.

Основная идея делать PPC - сравнить предсказанные значения с вашими данными, чтобы проверить, где они оба согласны, а где нет. Эта информация может использоваться, например, для улучшения модели. Выполнение

pm.sample_ppc(trace, model=linear_model, samples=100)

Дает вам 100 образцов, каждый из которых имеет 331 прогнозируемое наблюдение (так как в вашем примере y_tr имеет длину 331). Следовательно, вы можете сравнить каждую прогнозируемую точку данных с образцом размером 100, взятым из заднего. Вы получаете распределение прогнозируемых значений, потому что задний сам является распределением возможных параметров (распределение отражает неопределенность). Что касается аргументов sample_ppc: samples указать, сколько точек из заднего вы получаете, каждая точка является вектором параметров. size указывает, сколько раз вы используете этот вектор параметров для выборки прогнозируемых значений (по умолчанию size=1).

У вас есть больше примеров использования sample_ppc в этом учебнике