PyMC3 байесовской линейной регрессии предсказания с sklearn.datasets

Question

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

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

Я попробовал следующее: Generating predictions from inferred parameters in 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's 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 

Answer 1

Я думаю, что одна из проблем, с вашей моделью является то, что данные имеют очень разные масштабы, у вас есть ~ 0,3 диапазона для ваших «крестиков» и ~ 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); 

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

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 и вообще вы можете использовать те же прежде всего, для любых данных (что-то полезно, если у вас нет информационных ориентиров, которые вы можете использовать). Фактически, многие люди рекомендуют эту практику для обобщенных линейных моделей. Вы можете прочитать об этом в книге doing bayesian data analysis или Statistical Rethinking

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

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 в этом tutorial