La predizione bayesiana

La predizione bayesiana#

In questo capitolo, approfondiremo il concetto e l’importanza delle distribuzioni predittive a priori e a posteriori nel contesto dell’analisi di dataset. La distribuzione predittiva a posteriori è cruciale per verificare l’aderenza delle previsioni del modello ai dati osservati. Un allineamento tra queste previsioni e i dati effettivamente raccolti ci consente di convalidare l’accuratezza del modello nel rappresentare il processo generativo sottostante.

Parallelamente, la distribuzione predittiva a priori modella le aspettative dei dati prima di qualsiasi osservazione effettiva. Essa integra le nostre conoscenze preesistenti e le ipotesi sui parametri del modello, fornendo un quadro essenziale per la proiezione e l’interpretazione di fenomeni complessi in un contesto statistico avanzato. Questa distribuzione non solo predispone una struttura per l’analisi inferenziale, ma è anche fondamentale per la formulazione di nuove ipotesi e per la progettazione di esperimenti futuri.

Preparazione del Notebook#

%run ../_config/config.py # Import the configuration settings
import itertools
import logging
import statistics as stat
from scipy import stats
import cmdstanpy
from cmdstanpy import CmdStanModel
cmdstanpy.utils.get_logger().setLevel(logging.ERROR)
from scipy import stats
import seaborn as sns

La distribuzione predittiva a posteriori#

La distribuzione predittiva a posteriori (PPD) è un concetto essenziale nell’inferenza bayesiana e permette di valutare quanto bene un modello si adatti ai dati osservati. Secondo Gelman and Shalizi [GS13], questa metodologia fornisce una valutazione critica della coerenza tra i dati reali e quelli simulati dal modello. La verifica predittiva a posteriori (PPC) utilizza la PPD per confrontare direttamente i dati osservati con quelli generati dal modello, rivelando eventuali discrepanze che possono indicare problemi nella specificazione del modello. In pratica, la PPC funge da test diagnostico, consentendo di individuare e correggere eventuali lacune nel modello al fine di migliorarne le capacità predittive.

Per comprendere meglio il concetto, è utile considerare la distribuzione predittiva a posteriori in termini di un modello coniugato normale-normale. Supponiamo di voler predire la media di una distribuzione normale futura, basandoci sui dati osservati e sulle nostre conoscenze a priori. La PPD ci offre uno strumento per calcolare queste probabilità, combinando le informazioni provenienti dai dati osservati con quelle fornite dalla distribuzione a priori.

Ad esempio, immaginiamo di aver raccolto dati sulle altezze di 100 persone, ottenendo una media campionaria di 170 cm e una deviazione standard campionaria di 10 cm. Il nostro obiettivo è stimare la media delle altezze in un futuro campione di \(n=100\) persone. La nostra conoscenza a priori sulla media delle altezze è rappresentata da una distribuzione normale con media 175 cm e deviazione standard di 5 cm. In termini di notazione, possiamo esprimere questa distribuzione come \(P(\tilde{y}|\theta=\theta_1)\), dove \(\tilde{y}\) rappresenta un nuovo dato che è diverso dai dati attuali \(y\), e \(\theta_1\) è la media a posteriori. Tuttavia, in statistica bayesiana, è fondamentale incorporare tutta l’incertezza nei risultati. Poiché \(\theta_1\) è solo uno dei possibili valori per \(\theta\), dovremmo includere ogni valore di \(\theta\) per la nostra previsione. Per ottenere la migliore previsione, possiamo «mediare» le previsioni attraverso i diversi valori di \(\theta\), ponderando ciascun valore secondo la sua probabilità a posteriori.

La distribuzione risultante è la distribuzione predittiva a posteriori, che in notazione matematica è data da:

\[ P(\tilde{y}|y) = \int_\theta p(\tilde{y}|\theta, y) p(\theta|y) d\theta \]

In questo modo, la PPD combina le informazioni dai dati osservati con la conoscenza a priori, fornendo una previsione che riflette l’incertezza associata a tutti i possibili valori dei parametri del modello.

Distribuzione predittiva a posteriori nel modello normale-normale#

Nel modello coniugato normale-normale, se i dati osservati \(Y = \{y_1, y_2, ..., y_n\}\) sono modellati come provenienti da una distribuzione normale con media \(\mu\) e varianza \(\sigma^2\), e assumendo una distribuzione a priori normale per \(\mu\), la distribuzione a posteriori di \(\mu\) sarà anch’essa normale.

Formule della distribuzione predittiva a posteriori#

Dato che:

  1. I dati osservati \(y_i \sim \mathcal{N}(\mu, \sigma^2)\)

  2. La prior per \(\mu\) è \(\mu \sim \mathcal{N}(\mu_0, \tau_0^2)\)

La distribuzione a posteriori per \(\mu\) sarà:

\[ \mu \mid Y \sim \mathcal{N}(\mu_n, \tau_n^2) \]

dove:

\[ \mu_n = \frac{\tau_0^2 \bar{y} + \sigma^2 \mu_0}{\tau_0^2 + \sigma^2} \]

e

\[ \tau_n^2 = \frac{\tau_0^2 \sigma^2}{\tau_0^2 + \sigma^2} \]

Qui, \(\bar{y}\) è la media campionaria dei dati osservati.

Esempio pratico#

Consideriamo che:

  • \(\mu_0 = 175\) cm (media a priori)

  • \(\tau_0 = 5\) cm (deviazione standard a priori)

  • \(\bar{y} = 170\) cm (media campionaria)

  • \(\sigma = 10\) cm (deviazione standard campionaria)

  • \(n = 100\) (numero di osservazioni)

I parametri della distribuzione a posteriori sono:

\[ \mu_n = \frac{(5^2 \cdot 170) + (10^2 \cdot 175)}{5^2 + 10^2} = \frac{42500 + 175000}{25 + 100} = \frac{217500}{125} = 174 \quad \text{cm} \]
\[ \tau_n^2 = \frac{5^2 \cdot 10^2}{5^2 + 10^2} = \frac{2500}{125} = 20 \quad \text{cm}^2 \Rightarrow \tau_n = \sqrt{20} \approx 4.47 \quad \text{cm} \]

Pertanto, la distribuzione a posteriori per \(\mu\) è:

\[ \mu \mid Y \sim \mathcal{N}(174, 4.47^2) \]

Per la distribuzione predittiva a posteriori, dobbiamo considerare anche la varianza della distribuzione futura. Se stiamo predicendo per \(n_{\text{fut}}=100\) nuove osservazioni, la varianza della media predittiva sarà:

\[ \sigma_{\text{pred}}^2 = \tau_n^2 + \frac{\sigma^2}{n_{\text{fut}}} \]
\[ \sigma_{\text{pred}}^2 = 20 + \frac{10^2}{100} = 20 + 1 = 21 \quad \text{cm}^2 \Rightarrow \sigma_{\text{pred}} = \sqrt{21} \approx 4.58 \quad \text{cm} \]

Quindi, la distribuzione predittiva a posteriori è:

\[ \tilde{Y} \sim \mathcal{N}(174, 4.58^2) \]

Implementazione con cmdstanpy#

Per illustrare come viene generata la distribuzione predittiva a posteriori nel contesto del modello normale-normale, possiamo utilizzare cmdstanpy per eseguire l’analisi. Il codice seguente mostra come configurare il modello e generare previsioni.

# Dati osservati
y_observed = np.random.normal(170, 10, 100)
mean_y = np.mean(y_observed)
std_y = np.std(y_observed)

# Parametri a priori
mu_0 = 175
tau_0 = 5

# Parametri posteriori
tau_n_sq = (tau_0**2 * std_y**2) / (tau_0**2 + std_y**2)
tau_n = np.sqrt(tau_n_sq)
mu_n = (tau_0**2 * mean_y + std_y**2 * mu_0) / (tau_0**2 + std_y**2)

# Parametri predittivi
n_fut = 100
sigma_pred_sq = tau_n_sq + (std_y**2 / n_fut)
sigma_pred = np.sqrt(sigma_pred_sq)
mu_pred = mu_n

# Simulazioni
y_pred_samples = np.random.normal(mu_pred, sigma_pred, 1000)

# Grafico
plt.hist(y_pred_samples, bins=30, density=True, alpha=0.5, label='Posterior Predictive')
x = np.linspace(160, 190, 200)
plt.plot(x, stats.norm.pdf(x, mu_pred, sigma_pred), 'r-', lw=2, label='Predictive Distribution')
plt.axvline(x=mu_pred, color='k', linestyle='--', label='Mean Prediction')
plt.xlabel('Heights (cm)')
plt.ylabel('Density')
plt.title('Posterior Predictive Distribution for Heights')
plt.legend()
plt.show()
../_images/361a09e12b3680d802fa82e1a402d2c8836a2ecfe3d034819616969a67d48eef.png

Questo codice produce un grafico che illustra visivamente la distribuzione predittiva a posteriori per le altezze nel nostro campione di 100 nuove osservazioni, tenendo conto sia dei dati osservati che delle nostre aspettative iniziali.

In sintesi, la distribuzione predittiva a posteriori è stata generata nel modo seguente:

  1. Campioniamo un valore \(\mu\) dalla distribuzione a posteriori di \(\mu\).

  2. Campioniamo un valore \(\sigma\) dalla distribuzione a posteriori di \(\sigma\).

  3. Utilizziamo questi valori per generare un campione dalla distribuzione normale con parametri \(\mu\) e \(\sigma\).

  4. Ripetiamo questo processo molte volte.

La distribuzione dei valori ottenuti da questi campionamenti costituisce la distribuzione predittiva a posteriori.

Metodo MCMC#

Quando usiamo un PPL come Stan, la distribuzione predittiva viene stimata mediante il campionamento da una catena di Markov, che è particolarmente utile in scenari complessi dove l’analisi analitica potrebbe essere impraticabile. Attraverso i metodi MCMC, si stimano le potenziali osservazioni future \(p(\tilde{y} \mid y)\), indicate come \(p(y^{rep} \mid y)\), seguendo questi passaggi:

  1. Si campiona \(\theta_i \sim p(\theta \mid y)\): Viene selezionato casualmente un valore del parametro (o dei parametri) dalla distribuzione a posteriori.

  2. Si campiona \(y^{rep} \sim p(y^{rep} \mid \theta_i)\): Viene scelta casualmente un’osservazione dalla funzione di verosimiglianza, condizionata al valore del parametro (o dei parametri)ottenuto nel passo precedente.

Ripetendo questi due passaggi un numero sufficiente di volte, l’istogramma risultante approssimerà la distribuzione predittiva a posteriori.

Esaminiamo ora come ottenere la distribuzione predittiva a posteriori con Stan per i dati dell’esempio precedente. Iniziamo creando le distribuzioni a posteriori di \(\mu\) e \(\sigma\).

stan_ncp_file = os.path.join(
    project_directory, 'stan', 'gaussian_ncp.stan')

model_ncp = CmdStanModel(stan_file=stan_ncp_file)

print(model_ncp.code())

data {
    int<lower=1> N;
    vector[N] y;
}
transformed data {
    real y_mean = mean(y);
    real y_sd = sd(y);
}
parameters {
    real mu_raw;
    real<lower=0> sigma_raw;
}
transformed parameters {
    real mu;
    real<lower=0> sigma;
    mu = y_mean + y_sd * mu_raw;
    sigma = y_sd * sigma_raw;
}
model {
    // Priors:
    mu_raw ~ normal(0, 1);
    sigma_raw ~ normal(0, 1);
    // Likelihood:
    y ~ normal(mu, sigma);
}
generated quantities {
    vector[N] y_rep;
    for (n in 1:N) {
        y_rep[n] = normal_rng(mu, sigma);
    }
}

Definiamo un dizionario che contiene i dati.

stan_data = {
    'N': len(y_observed), 
    'y': y_observed
}

print(stan_data)

{'N': 100, 'y': array([171.96581381, 164.47710641, 157.91426703, 158.32818828,
       163.86028418, 173.36934944, 154.07963372, 164.40764344,
       163.18301734, 174.07557036, 179.96191717, 189.1043359 ,
       172.31037122, 166.12502256, 160.60694872, 171.15798161,
       158.69657965, 157.21990299, 178.62771666, 160.35962282,
       162.28151759, 163.61569088, 151.0949321 , 169.40645352,
       147.87575638, 166.19409576, 149.31394323, 154.11754244,
       152.62083154, 182.28170464, 159.26720722, 171.56296352,
       182.05327316, 166.91587353, 175.63514612, 184.93380777,
       172.36334656, 169.90031078, 176.81658641, 163.5765079 ,
       158.87766777, 157.39365575, 172.60333912, 152.1230542 ,
       167.80733883, 165.8121181 , 173.33916904, 185.46442278,
       161.11319673, 178.7186444 , 158.86008323, 158.94742562,
       156.95008148, 157.82271689, 150.84899944, 163.98863641,
       184.39146361, 174.82034233, 171.60387418, 163.36719089,
       175.18099035, 175.86158476, 174.92277822, 190.25430142,
       160.47608042, 173.25373406, 165.23459932, 180.50866441,
       165.4129635 , 187.15268685, 171.73568326, 169.71540944,
       184.28396189, 172.25606125, 162.44021921, 163.83979063,
       172.20506434, 162.36493187, 170.83325995, 182.17875676,
       165.60317875, 170.80602341, 175.32058725, 181.30082987,
       165.60142294, 159.39040688, 179.14782863, 198.8319763 ,
       178.3322677 , 160.05915761, 179.38638789, 180.81940259,
       162.99635058, 169.45240358, 165.01596354, 174.44816503,
       175.2240379 , 179.48196389, 157.49157985, 165.47361839])}

Eseguiamo il campionamento MCMC:

trace_ncp = model_ncp.sample(
    data=stan_data,
    iter_warmup = 1_000,
    iter_sampling = 2_000,
    seed = 123,
    show_progress = False, 
    show_console = False
)

Un sommario delle distribuzioni a posteriori dei parametri si ottiene nel modo seguente:

az.summary(trace_ncp, var_names=['mu', 'sigma'], round_to=2)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu 168.82 1.01 166.97 170.74 0.01 0.01 6613.96 5236.04 1.0
sigma 10.20 0.72 8.86 11.53 0.01 0.01 7101.57 5186.43 1.0

Convertiamo l’oggetto sample_ncp in un oggetto di classe InferenceData:

idata = az.from_cmdstanpy(
    posterior=trace_ncp, 
    posterior_predictive='y_rep', 
    observed_data={"y": y_observed}
)

La distribuzione predittiva a posteriori è utilizzata per eseguire i controlli predittivi a posteriori (PPC), noti come Posterior Predictive Checks. I PPC consistono in un confronto grafico tra \(p(y^{rep} \mid y)\), ossia la distribuzione delle osservazioni future previste, e i dati osservati \(y\). Questo confronto visivo permette di valutare se il modello utilizzato è adeguato per descrivere le proprietà dei dati osservati.

Oltre al confronto grafico tra le distribuzioni \(p(y)\) e \(p(y^{rep})\), è possibile effettuare un confronto tra le distribuzioni di varie statistiche descrittive calcolate su diversi campioni \(y^{rep}\) e le corrispondenti statistiche calcolate sui dati osservati. Tipicamente, vengono considerate statistiche descrittive come la media, la varianza, la deviazione standard, il minimo o il massimo, ma è possibile confrontare qualsiasi altra statistica rilevante.

I controlli predittivi a posteriori offrono un valido strumento per un’analisi critica delle prestazioni del modello e, se necessario, per apportare eventuali modifiche o considerare modelli alternativi più adatti ai dati in esame.

Possiamo ora usare ArviZ per generare il posterior-predictive plot:

_ = az.plot_ppc(idata, data_pairs={"y": "y_rep"}, num_pp_samples=500)
../_images/2aef173908819bf27e6a9a1434df9926717469de5361ae9906837880cedf6ec0.png

Distribuzione Predittiva a Priori#

Le verifiche predittive a priori generano dati utilizzando unicamente le distribuzioni a priori, ignorando i dati osservati, al fine di valutare se tali distribuzioni a priori sono appropriate (Gabry et al. 2019). La distribuzione predittiva a priori è quindi simile alla distribuzione predittiva a posteriori, ma senza dati osservati, rappresentando il caso limite di una verifica predittiva a posteriori senza dati.

Questo processo può essere realizzato facilmente simulando i parametri secondo le distribuzioni a priori e poi simulando i dati in base al modello dati i parametri simulati. Il risultato è una simulazione dalla distribuzione congiunta, che è quindi una simulazione dalla distribuzione predittiva a priori.

Questa procedura è fondamentale per verificare se le ipotesi a priori sono realistiche e adeguate prima di raccogliere o utilizzare i dati osservati. Se i dati simulati dalla distribuzione predittiva a priori non risultano plausibili, potrebbe essere necessario rivedere le scelte delle distribuzioni a priori.

Consideriamo, quale esempio, il caso discusso in precedenza di un campione di dati gaussiani e di un modello gaussiano in cui le distribuzioni a priori per μ e σ sono gaussiane.

stan_file = os.path.join(
    project_directory, 'stan', 'gaussian_model_prior.stan')

model_gauss = CmdStanModel(stan_file=stan_file)

print(model_gauss.code())

data {
  int<lower=0> N;            // number of observations
}
generated quantities {
  real mu;                   // parameter of interest
  real<lower=0> sigma;       // known standard deviation of y
  array[N] real y_rep;       // prior predictive samples

  // Priors
  mu = normal_rng(175, 5);
  sigma = 10;

  // Generate prior predictive samples
  for (n in 1:N) {
    y_rep[n] = normal_rng(mu, sigma);
  }
}

Il codice precedente illustra come definire il modello Stan per generare campioni predittivi a priori. In questo esempio, mu e sigma sono generati dalle loro rispettive distribuzioni a priori e usati per generare campioni di dati simulati y_rep.

# Stan data dictionary
stan_data = {
    'N': len(y_observed)
}

prior_predictive_samples = model_gauss.sample(
    data=stan_data, 
    fixed_param=True, 
    iter_sampling=1000, 
    iter_warmup=1, 
    chains=1,
    show_progress=False, 
    show_console=False
)

# Extract the relevant variables
y_rep_samples = prior_predictive_samples.stan_variable('y_rep')
y_rep_flattened = y_rep_samples.flatten()

# Check the statistics of y_rep values
y_rep_mean = np.mean(y_rep_flattened)
y_rep_std = np.std(y_rep_flattened)

print(f'Mean of y_rep: {y_rep_mean}')
print(f'Standard Deviation of y_rep: {y_rep_std}')

Mean of y_rep: 175.04561795
Standard Deviation of y_rep: 11.105448274627541

# Create a KDE plot
plt.figure(figsize=(10, 6))
sns.kdeplot(y_rep_flattened, fill=True)
plt.title('Prior Predictive Check')
plt.xlabel('Height (cm)')
plt.ylabel('Density')
plt.axvline(x=y_rep_mean, color='r', linestyle='--', label=f'Mean: {y_rep_mean:.2f}')
plt.legend()
plt.show()
../_images/7436dbe6f506fb581996d8487f6cc2fc750dc1639229fc86edffd6404fb34248.png

Questo approccio assicura che le distribuzioni a priori siano realistiche e adeguate, permettendo di identificare e correggere eventuali ipotesi errate prima di procedere con l’analisi dei dati osservati.

Considerazioni conclusive#

Le distribuzioni predittive a priori e a posteriori sono generate in maniera simile, con la seguente differenza.

  • Distribuzione Predittiva a Priori: Questa distribuzione rappresenta le nostre previsioni sui dati prima di osservare qualsiasi dato effettivo. In questo caso, prendiamo valori dei parametri dalla distribuzione a priori e li utilizziamo nella funzione di verosimiglianza per generare dati predittivi. La distribuzione di questi dati generati è la nostra distribuzione predittiva a priori, che riflette le nostre conoscenze e incertezze prima dell’osservazione dei dati.

  • Distribuzione Predittiva a Posteriori: Dopo aver osservato i dati, aggiorniamo le nostre credenze sulla distribuzione dei parametri usando il teorema di Bayes, ottenendo così la distribuzione a posteriori dei parametri. La distribuzione predittiva a posteriori viene generata prendendo valori dei parametri dalla distribuzione a posteriori (che incorpora le informazioni dai dati osservati) e inserendoli nella funzione di verosimiglianza per generare nuovi dati predittivi. Questa distribuzione riflette le nostre previsioni sui dati futuri o non osservati, dopo aver considerato i dati attuali.

La differenza chiave tra le due distribuzioni predittive è quindi la distribuzione dei parametri utilizzata per generare i dati: il prior nel caso della distribuzione predittiva a priori, e il posterior nel caso della distribuzione predittiva a posteriori. La distribuzione predittiva a posteriori è generalmente più informativa perché tiene conto dei dati osservati.

È fondamentale, per l’integrità del modello, che la distribuzione predittiva a posteriori rifletta la distribuzione dei dati osservati. Per validare questa corrispondenza, si utilizzano le verifiche predittive a posteriori, confrontando la distribuzione predittiva con i dati empirici tramite stime di densità Kernel (KDE). Questo confronto consente di valutare l’efficacia del modello nell’approssimare la struttura sottostante dei dati e la sua capacità di guidare previsioni affidabili.

Informazioni sull’Ambiente di Sviluppo#

%load_ext watermark
%watermark -n -u -v -iv -w -m -p jax

Last updated: Sun Jun 16 2024

Python implementation: CPython
Python version       : 3.12.3
IPython version      : 8.25.0

jax: 0.4.27

Compiler    : Clang 16.0.6 
OS          : Darwin
Release     : 23.4.0
Machine     : arm64
Processor   : arm
CPU cores   : 8
Architecture: 64bit

cmdstanpy : 1.2.3
arviz     : 0.18.0
seaborn   : 0.13.2
pandas    : 2.2.2
numpy     : 1.26.4
matplotlib: 3.8.4
logging   : 0.5.1.2
scipy     : 1.13.1

Watermark: 2.4.3
Contenuti