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
.PyMC3 bayesiana regressione lineare previsione con sklearn.datasets
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
suona bene, io sicuramente capito. lo porterò via adesso –
Fatto già, e grazie! – halfer