PyMC3 bayesiana regressione lineare previsione con sklearn.datasets

Question

ho cercato di implementare bayesiano di regressione lineare modelli utilizzando PyMC3 con REALE DATI (vale a dire non dalla funzione lineare + rumore gaussiano) dai set di dati in sklearn.datasets. Ho scelto il set di dati di regressione con il numero minimo di attributi (ad esempio load_diabetes()) la cui forma è (442, 10); ovvero, 442 samples e 10 attributes.

Credo di avere il modello funzionante, i posteriori sembrano abbastanza decenti da cercare di prevedere per capire come funziona questa roba ma ... Mi sono reso conto che non ho idea di come prevedere con questi modelli bayesiani! Sto cercando di evitare di utilizzare la notazione glm e patsy perché è difficile per me capire cosa sta effettivamente accadendo quando si utilizza quello.

ho cercato seguente: Generating predictions from inferred parameters in pymc3 e anche http://pymc-devs.github.io/pymc3/posterior_predictive/ ma il mio modello è o estremamente terribile a prevedere o che sto facendo male.

Se in realtà sto facendo correttamente la previsione (che probabilmente non sono), allora qualcuno può aiutarmi a ottimizzare il mio modello. Non so se almeno mean squared error, absolute error, o qualcosa del genere funziona in strutture bayesiane. Idealmente, mi piacerebbe ottenere un array di number_of_rows = la quantità di righe nel mio set di prova attributo/dati X_te e il numero di colonne da campionare dalla distribuzione posteriore.

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

Credo che uno dei problemi con il vostro modello è che i dati ha scale molto diverse, avete ~ 0,3 gamma per il vostro "xs" e ~ 300 per il vostro "Ys ". Quindi dovresti aspettarti pendenze maggiori (e sigma) che i tuoi priori stanno specificando. Un'opzione logica è quella di regolare i priori, come nell'esempio seguente.

#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); 

posteriori controlli predittivi, dimostra che hai un modello più o meno ragionevole.

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

Un'altra opzione sarà quella di inserire i dati nella stessa scala standardizzando esso, così facendo si otterrà che la pendenza dovrebbe essere di circa + -1 e, in generale, è possibile utilizzare lo stesso anteriore diffuso per qualsiasi dato (qualcosa di utile a meno che tu non abbia dei priori informativi che puoi usare). In effetti, molte persone raccomandano questa pratica per i modelli lineari generalizzati. Il tuo grado di saperne di più su questo nel libro doing bayesian data analysis o Statistical Rethinking

Se si vuole prevedere i valori Sono disponibili diverse opzioni di uno è quello di utilizzare la media dei parametri di dedurre, come:

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

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

Un'altra opzione è quella di utilizzare pm.sample_ppc per ottenere campioni di valori previsti che tengano conto dell'incertezza delle stime.

L'idea principale di fare PPC è confrontare i valori previsti con i dati per verificare dove sono entrambi d'accordo e dove no. Questa informazione può essere utilizzata ad esempio per migliorare il modello. Facendo

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

vi darà 100 campioni ciascuna con 331 predetto osservazione (in quanto nel tuo esempio y_tr ha lunghezza 331). Quindi è possibile confrontare ciascun punto dati previsto con un campione di dimensione 100 prelevato dal posteriore. Si ottiene una distribuzione dei valori previsti perché il posteriore è di per sé una distribuzione di parametri possibili (la distribuzione riflette l'incertezza). Per quanto riguarda gli argomenti di sample_ppc: samples specificare quanti punti dalla parte posteriore si ottiene, ogni punto è vettore di parametri. size specifica quante volte si utilizza quel vettore di parametri per campionare i valori previsti (per impostazione predefinita size=1).

avete più esempi di utilizzo sample_ppc in questo tutorial