przewidywania PyMC3 Bayesa regresji liniowej z sklearn.datasets

Question

Próbuję zaimplementować Bayesa regresji liniowej modele używając PyMC3 z rzeczywistych danych (to znaczy nie od funkcji liniowej + szum Gaussa) z zestawów danych w sklearn.datasets. Wybrałem zbiór danych regresji o najmniejszej liczbie atrybutów (tj. load_diabetes()), którego kształt to (442, 10); to jest 442 samples i .

Wierzę, że mam działający model, a posteriors wyglądają na dość przyzwoicie, aby spróbować przewidzieć, jak to działa, ale ... Zdałem sobie sprawę, że nie mam pojęcia, jak przewidzieć z tymi modelami Bayesian! Staram się unikać notacji glm i patsy, ponieważ trudno mi zrozumieć, co faktycznie dzieje się podczas korzystania z tego.

Próbowałem następujące: Generating predictions from inferred parameters in pymc3 a także http://pymc-devs.github.io/pymc3/posterior_predictive/ ale mój model jest albo bardzo straszny w przewidywaniu czy robię to źle.

Jeśli właściwie wykonuję przewidywanie (co prawdopodobnie nie jest), to każdy może mi pomóc zoptymalizować mój model. Nie wiem, czy najmniej mean squared error, absolute error, lub coś podobnego działa w frameworkach Bayesian. Idealnie chciałbym uzyskać tablicę number_of_rows = ilość wierszy w moim zestawie testów atrybutów/danych X_te oraz liczbę kolumn, które mają być próbkami z rozkładu wstecznego.

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

Myślę, że jednym z problemów z modelu jest to, że dane mają bardzo różne skale, masz ~ 0,3 zakres dla „xs” i ~ 300 dla „Ys ". Dlatego powinieneś oczekiwać większych stoków (i sigmy), które określają twoi priors. Jedną z logicznych opcji jest dostosowanie twoich wartości początkowych, jak w poniższym przykładzie.

#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 sprawdza predykcyjne, pokazuje, że masz więcej lub mniej rozsądny model.

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') 

Inną opcją będzie umieścić dane w tej samej skali przez standaryzację go, robiąc tak dostaniesz, że nachylenie powinno wynosić około + -1 iw ogóle można użyć tego samego rozproszony wcześniej dla jakichkolwiek danych (coś pożytecznego, chyba że masz informacje, których możesz użyć). W rzeczywistości wiele osób zaleca tę praktykę w przypadku uogólnionych modeli liniowych. Twój może przeczytać więcej na ten temat w książce doing bayesian data analysis lub Statistical Rethinking

Jeśli chcesz przewidzieć wartości masz kilka opcji, jedna jest użycie średnia z uzyskanych danych parametrów, takich jak:

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

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

Innym rozwiązaniem jest użyj pm.sample_ppc, aby pobrać próbki przewidywanych wartości, które uwzględniają niepewność oszacowań.

Główną ideą wykonywania PPC jest porównanie przewidywanych wartości z danymi w celu sprawdzenia, gdzie oboje się zgadzają, a gdzie nie. Informacje te można wykorzystać na przykład do ulepszenia modelu. Robi

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

daje 100 próbek każdy z 331 przewidział obserwację (ponieważ w swoim przykładzie y_tr ma długość 331). W związku z tym można porównać każdy przewidywany punkt danych z próbką o wielkości 100 pobraną z tyłu. Otrzymujesz rozkład przewidywanych wartości, ponieważ sam tylny jest rozkładem możliwych parametrów (rozkład odzwierciedla niepewność). Odnośnie do argumentów z sample_ppc: samples określ, ile punktów z tyłu otrzymujesz, każdy punkt jest wektorem parametrów. size określa, ile razy użyjesz wektora parametrów do wypróbowania przewidywanych wartości (domyślnie size=1).

Masz więcej przykładów wykorzystania sample_ppc w tym tutorial