Inferenza bayesiana con PyMC

Inferenza bayesiana con PyMC#

In questo capitolo, faremo uso di PyMC, una libreria Python specializzata nella modellazione statistica bayesiana. A differenza del capitolo precedente, in cui abbiamo utilizzato le funzioni di scipy.stats, qui ricorreremo a un Linguaggio di Programmazione Probabilistico (PPL) per implementare il nostro campionatore. L’uso di PPL semplifica notevolmente l’attuazione dell’inferenza bayesiana, permettendo agli utenti di focalizzarsi sulla costruzione del modello e liberandoli da dettagli matematici e computazionali.

Nel contesto di questa metodologia avanzata, in questo e nei successivi capitoli esploreremo i modelli a un parametro. Questi modelli rappresentano una classe di distribuzioni campionarie identificate da un unico parametro sconosciuto. In particolare, ci concentreremo sull’inferenza bayesiana applicata a due tipi di modelli a un parametro: il modello binomiale e il modello di Poisson. Oltre ad essere utili strumenti statistici, questi modelli forniscono anche un ambiente semplice all’interno del quale è possibile apprendere le basi dell’analisi dei dati bayesiana, incluse le distribuzioni a priori coniugate, le distribuzioni predittive e le regioni di confidenza.

Campionatore No-U-Turn#

L’algoritmo Metropolis, che abbiamo esaminato in precedenza, offre un metodo per generare campioni da distribuzioni di probabilità attraverso la creazione di una catena di Markov, con la distribuzione desiderata come distribuzione di equilibrio (o stazionaria). Tuttavia, per modelli complessi, tale algoritmo può risultare inefficiente e richiedere un tempo considerevole per convergere verso una distribuzione stazionaria. Per affrontare questa sfida, sono stati sviluppati algoritmi Monte Carlo a catena di Markov (MCMC) più efficienti, tra cui il campionatore No-U-Turn (NUTS), gli algoritmi Metropolis-Hastings, il campionatore Gibbs e il campionatore Hamiltonian Monte Carlo. Questi algoritmi sono stati implementati in vari PPL, tra cui PyMC e Stan [CGH+17].

Nel presente capitolo useremo PyMC per condurre l’inferenza sulle proporzioni binomiali.

Preparazione del Notebook#

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pymc as pm
import arviz as az
import xarray as xr
import scipy.stats as stats
import warnings

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=Warning)

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 42
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")
sns.set_theme(palette="colorblind")

Un esempio pratico di utilizzo del linguaggio probabilistico PyMC consiste nell’eseguire una semplice operazione aritmetica. Iniziamo a sommare due numeri interi in Python.

# adding 2 integers in Python
a = 2
b = 3
c = a + b
print(c)

Facciamo ora la stessa cosa usando PyMC.

# adding 2 random variables in PyMC
with pm.Model() as example:
    a = pm.Normal("a", 2, 0.5)
    b = pm.Normal("b", 3, 0.2)
    c = pm.Deterministic("c", a + b)
    trace_1 = pm.sample()

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [a, b]

100.00% [8000/8000 00:09<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 10 seconds.

with pm.Model() as example:

In questa linea, viene creato un nuovo modello probabilistico PyMC utilizzando un contesto Python with. All’interno di questo contesto, tutte le variabili casuali e le operazioni definite saranno automaticamente associate al modello chiamato example.

a = pm.Normal("a", 2, 0.5)

Qui viene definita una variabile casuale normalmente distribuita con un valore medio (\(\mu\)) di 2 e una deviazione standard (\(\sigma\)) di 0.5. Questa variabile è denominata “a”.

b = pm.Normal("b", 3, 0.2)

Qui viene definita un’altra variabile aleatoria normalmente distribuita, ma questa volta con un valore medio di 3 e una deviazione standard di 0.2. Questa variabile è denominata “b”.

c = pm.Deterministic("c", a + b)

In questa linea, viene creata una nuova variabile deterministica “c” che è definita come la somma delle variabili “a” e “b”. Essendo una variabile deterministica, il suo valore è completamente determinato dai valori di “a” e “b”.

trace_1 = pm.sample()

Infine, questa linea avvia il processo di campionamento MCMC (Monte Carlo Markov Chain) per generare un campione dalla distribuzione a posteriori del modello. Il risultato del campionamento è salvato in un oggetto chiamato trace_1, che contiene le tracce delle variabili aleatorie durante il campionamento.

Per rappresentare graficamente i risultati del processo di campionamento, possiamo utilizzare la funzione plot_trace() fornita dalla libreria ArviZ. Questa funzione crea un diagramma delle tracce che visualizza sia la serie temporale dei valori campionati (traccia) sia la loro distribuzione.

az.plot_trace(trace_1, figsize=(10, 6), combined=True)
plt.tight_layout()

../_images/f526fd8bb83bc8822d83b39f1d005b10cb012072393dfccab753cb1adb329680.png

Il grafico delle tracce illustra che la moda della distribuzione a posteriori della variabile \(c\) corrisponde al risultato che prevedevamo.

Generazione X#

Ora analizzeremo un caso di inferenza basata su dati reali, in contrapposizione alla situazione deterministica esaminata precedentemente. Prenderemo nuovamente in esame i dati discussi nel capitolo precedente, relativi agli artisti della Generazione X rappresentati al MOMA. Ricordiamo che abbiamo registrato 14 successi su 100 tentativi e abbiamo basato il parametro \( \theta \) (la probabilità di appartenere alla Generazione X o a generazioni successive) su una distribuzione Beta(4, 6).

Il nostro modello, quindi, si articola come segue:

\[\begin{split} \begin{align*} \theta &\sim \text{Beta}(4, 6), \\ y &\sim \text{Bin}(n, \theta). \end{align*} \end{split}\]

In questa impostazione, la prima riga definisce la distribuzione a priori, \(p(\theta)\), mentre la seconda riga determina la funzione di verosimiglianza, \(p(y_1, \ldots, y_{100} | \theta)\).

Nel caso presente abbiamo deciso di rappresentare la nostra incertezza a priori rispetto al parametro \( \theta \) mediante una Beta(4, 6).

x = np.linspace(0, 1, 1000)

prior_density = stats.beta.pdf(x, 4, 6)

plt.fill_between(x, prior_density, alpha=0.5)
plt.xlabel("Parameter Value")
plt.ylabel("Density")
plt.title("Prior Density");

../_images/b7967d6e8a53ff9c0ca5136450c3d7ffdeb960f397b230886b3064cbe02790db.png

Per questa distribuzione a priori e il modello di campionamento sopra indicato, la regola di Bayes fornisce

\[\begin{split} \begin{align*} p(\theta \mid y_1, \ldots, y_{100}) &= \frac{p(y_1, \ldots, y_{100} \mid \theta) p(\theta)}{p(y_1, \ldots, y_{100})} \\ &\propto p(y_1, \ldots, y_{100} \mid \theta) p(\theta) \\ &\propto \prod_{i=1}^{100} p(y_i \mid \theta) p(\theta), \end{align*} \end{split}\]

dove \(p(y_i \mid \theta)\) è la distribuzione bernoulliana

\[ p(Y = y \mid \theta) = \theta^y (1-\theta)^{1-y} \]

e \(p(\theta)\) è una distribuzione Beta(4,6).

Nel contesto specifico che stiamo analizzando, è possibile trovare una soluzione analitica per stimare la distribuzione a posteriori. Per esempio, se \( y = 14 \) e \( n = 100 \), la distribuzione a posteriori sarà una Beta(18, 92). Infatti, la distribuzione Beta è un “prior coniugato” per la distribuzione binomiale, il che significa che se iniziamo con un prior Beta e osserviamo dati che seguono una distribuzione binomiale, il posterior sarà ancora una distribuzione Beta.

La formula per aggiornare i parametri della distribuzione Beta a posteriori \( \text{Beta}(\alpha', \beta') \) a partire da un prior \( \text{Beta}(\alpha, \beta) \) e dati binomiali con \( y \) successi e \( n \) prove è la seguente:

\[ \alpha' = \alpha + y \]

\[ \beta' = \beta + n - y \]

Nel nostro caso, \( \alpha = 4 \), \( \beta = 6 \), \( y = 14 \), e \( n = 100 \). Utilizzando queste formule:

\[ \alpha' = 4 + 14 = 18 \]

\[ \beta' = 6 + 100 - 14 = 92 \]

otteniamo che la distribuzione a posteriori è una \( \text{Beta}(18, 92) \).

Sebbene in questo caso particolare possiamo facilmente ottenere una soluzione analitica per la distribuzione a posteriori, tale opportunità è rara nella maggior parte dei modelli di inferenza bayesiana. Spesso, per gestire modelli più complessi, è necessario ricorrere a tecniche di approssimazione numerica, come i metodi MCMC. Nell’esempio qui discusso, applicheremo i metodi MCMC a un caso per il quale abbiamo già una soluzione analitica, permettendoci così di confrontare i risultati derivati dall’approccio analitico con quelli ottenuti attraverso l’approssimazione numerica.

Dedurre una proporzione con PyMC#

Adesso eseguiremo l’analisi che abbiamo precedentemente condotto utilizzando l’algoritmo di Metropolis usando PyMC. Assumiamo di avere già installato il pacchetto PyMC. Una volta installato, dovremo importare le librerie necessarie, tra cui Matplotlib, Numpy, Scipy, Arviz e ovviamente PyMC stesso.

Esamineremo ora come specificare il modello beta-binomiale attraverso PyMC. Per condurre l’analisi tramite PyMC, sarà prima necessario delineare la struttura del modello bayesiano e, successivamente, eseguire il campionamento dalla distribuzione a posteriori. Approfondiremo entrambi questi passaggi nell’ambito del nostro esempio.

Dati#

y = 14
ntrials = 100

Distribuzione a priori#

I parametri della distribuzione Beta, scelta come distribuzione a priori per \(\theta\), sono i seguenti.

alpha_prior = 4
beta_prior = 6

Panoramica dell’architettura di PyMC#

Moduli e componenti#

PyMC è organizzato in una serie di moduli che, pur essendo quasi indipendenti, sono interconnessi e supportano la modellazione bayesiana. Questi moduli includono componenti fondamentali come modelli e distribuzioni di probabilità, vari algoritmi di inferenza e anche elementi più avanzati.

Distribuzioni#

PyMC dispone di un’ampia gamma di distribuzioni di probabilità, sia quelle comunemente usate che altre più specializzate, tutte implementate come classi separate.

Gestione del modello#

Il modulo Model offre la classe Model, che ingloba tutti gli aspetti del modello bayesiano definito dall’utente. L’interazione con un’istanza del modello si svolge tramite un gestore di contesto, che facilita l’aggiunta automatica di variabili al modello e l’instaurazione di relazioni tra esse, come le relazioni genitore-figlio tra variabili aleatorie.

Step Methods#

Per le distribuzioni continue, il campionatore predefinito è basato sul No-U-Turn Sampler (NUTS). Tuttavia, gli utenti hanno la possibilità di sovrascrivere manualmente questo algoritmo predefinito.

Campionamento e funzioni ausiliarie#

PyMC non solo supporta vari metodi MCMC, ma permette anche il campionamento dalle distribuzioni predittive a priori e a posteriori. Inoltre, offre la flessibilità di utilizzare la stessa definizione del modello sia per calcolare le distribuzioni a posteriori che per quelle predittive, senza necessità di ulteriori interventi da parte dell’utente.

Commenti e considerazioni finali#

Questo capitolo ha illustrato come utilizzare PyMC per ottenere la distribuzione a posteriori quando ci troviamo di fronte a un caso di distribuzione beta-binomiale, caratterizzato dalle relazioni di verosimiglianza binomiale e dalla distribuzione a priori Beta. Inoltre, è stata condotta un’analisi comparativa tra la soluzione ottenuta tramite PyMC e la soluzione analitica del problema discusso, dimostrando la convergenza tra i due approcci.

%load_ext watermark
%watermark -n -u -v -iv -w

Last updated: Tue Jan 23 2024

Python implementation: CPython
Python version       : 3.11.7
IPython version      : 8.19.0

xarray    : 2023.12.0
pymc      : 5.10.3
scipy     : 1.11.4
matplotlib: 3.8.2
seaborn   : 0.13.0
numpy     : 1.26.2
arviz     : 0.17.0

Watermark: 2.4.3

Inferenza bayesiana con PyMC

Contents

Inferenza bayesiana con PyMC#

Campionatore No-U-Turn#

Preparazione del Notebook#

Generazione X#

Dedurre una proporzione con PyMC#

Dati#

Distribuzione a priori#

Creazione del Modello#

Inizializzazione#

Variabili aleatorie e distribuzioni#

Variabili osservate e non osservate#

Relazioni gerarchiche#

Contenitori di dati#

Inferenza#

Metodi di inferenza in PyMC: da Metropolis a NUTS#

Utilizzo dell’Algoritmo Metropolis#

L’oggetto `InferenceData` in PyMC#

Caratteristiche di un `xarray.Dataset`#

Panoramica dell’architettura di PyMC#

Moduli e componenti#

Distribuzioni#

Gestione del modello#

Step Methods#

Campionamento e funzioni ausiliarie#

Commenti e considerazioni finali#

Inferenza bayesiana con PyMC

Contents

Inferenza bayesiana con PyMC#

Campionatore No-U-Turn#

Preparazione del Notebook#

Generazione X#

Dedurre una proporzione con PyMC#

Dati#

Distribuzione a priori#

Creazione del Modello#

Inizializzazione#

Variabili aleatorie e distribuzioni#

Variabili osservate e non osservate#

Relazioni gerarchiche#

Contenitori di dati#

Inferenza#

Metodi di inferenza in PyMC: da Metropolis a NUTS#

Utilizzo dell’Algoritmo Metropolis#

L’oggetto InferenceData in PyMC#

Caratteristiche di un xarray.Dataset#

Panoramica dell’architettura di PyMC#

Moduli e componenti#

Distribuzioni#

Gestione del modello#

Step Methods#

Campionamento e funzioni ausiliarie#

Commenti e considerazioni finali#

L’oggetto `InferenceData` in PyMC#

Caratteristiche di un `xarray.Dataset`#