Monte Carlo a Catena di Markov

In precedenza, abbiamo esplorato diversi esempi di inferenza bayesiana relativi alla distribuzione a posteriori di un singolo parametro, come nel caso del modello bernoulliano. Abbiamo anche discusso l’utilizzo di approcci come l’approssimazione tramite griglia e i metodi dei priori coniugati per ottenere o approssimare la distribuzione a posteriori. In questo capitolo, ci concentreremo sul metodo di simulazione e spiegheremo perché sono necessari approcci speciali noti come metodi di Monte Carlo a Catena di Markov (MCMC).

Preparazione del Notebook

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy as sp
import scipy.stats as stats
import statsmodels.api as sm
from statsmodels.graphics import tsaplots
import arviz as az
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")

Il denominatore bayesiano

Nel campo dell’approccio bayesiano, il nostro obiettivo primario si concentra sulla determinazione della distribuzione a posteriori \(p(\theta \mid y)\) di un parametro \(\theta\) (nel caso più semplice), utilizzando come base sia i dati osservati \(y\) che la distribuzione a priori \(p(\theta)\). Questo processo si avvale del teorema di Bayes per calcolare tale distribuzione.

\[ p(\theta \mid y) = \frac{p(y \mid \theta) p(\theta)}{p(y)}. \]

Tuttavia, ci imbattiamo spesso in una sfida significativa: il calcolo dell’evidenza \(p(y)\) può rivelarsi estremamente complesso, specialmente per modelli più articolati, rendendo difficile ottenere con precisione la distribuzione a posteriori.

Una soluzione possibile si trova nelle distribuzioni a priori coniugate, che offrono un metodo analitico per determinare la distribuzione a posteriori. Questo però limita la selezione delle distribuzioni a priori e di verosimiglianza. Per superare questa limitazione, soprattutto in modelli dove i metodi di campionamento a griglia non sono applicabili, si utilizzano i Metodi di Catena di Markov a Monte Carlo (MCMC). Questi metodi rappresentano una potente alternativa, consentendo di derivare la distribuzione a posteriori basandosi su presupposti teorici e senza restrizioni nella scelta delle distribuzioni.

L’approccio Monte Carlo si basa sulla generazione di sequenze di numeri casuali per creare un campione ampio di osservazioni proveniente dalla distribuzione a posteriori. Da questi campioni, possiamo poi stimare empiricamente le proprietà di interesse. Questo approccio richiede l’uso di metodi computazionalmente intensivi, e con la crescente potenza di calcolo dei computer moderni, tali metodi stanno diventando sempre più accessibili e popolari nell’analisi dei dati.

Catene di Markov Monte Carlo

Nei capitoli precedenti abbiamo già esplorato l’efficacia della simulazione nel campo della teoria delle probabilità. Un esempio classico è il problema di Monty Hall, che mostra come la simulazione ripetuta possa fornire stime affidabili per la media e la varianza di variabili casuali. Inoltre, applicando la legge dei grandi numeri, queste stime diventano più accurate all’aumentare del numero di simulazioni. Ciò evidenzia la potenza dei metodi Monte Carlo, che evitano la necessità di calcolare somme o integrali complessi.

Il Metodo di Monte Carlo

Il Metodo di Monte Carlo, sviluppato durante il Progetto Manhattan negli anni ‘40, utilizza numeri casuali per risolvere problemi matematici complessi. Originariamente ideato da Stanislaw Ulam e successivamente implementato da John von Neumann, il metodo prende il nome dallo zio giocatore d’azzardo di Ulam, come suggerito da Nicholas Metropolis. Da allora, questo metodo è diventato una tecnica fondamentale in diverse discipline, contribuendo significativamente alla risoluzione di problemi complessi.

La metodologia di Monte Carlo genera un’ampia serie di punti casuali per stimare quantità di interesse, come l’integrazione numerica. Un esempio classico è l’approssimazione dell’integrale di un cerchio in 2D, dove il rapporto tra il numero di punti che cadono all’interno del cerchio e tutti i campioni fornisce un’approssimazione dell’area (per un esempio numerico, si veda Simulazione Monte Carlo).

Per illustrare ulteriormente questo concetto, consideriamo ora una distribuzione continua \(p(\theta \mid y)\) con una media \(\mu\). Se siamo in grado di generare una sequenza di campioni casuali \(\theta^{(1)}, \theta^{(2)}, \dots, \theta^{(T)}\) indipendenti e identicamente distribuiti secondo \(p(\theta \mid y)\), possiamo stimare il valore atteso teorico di \(\theta\) utilizzando la media campionaria \(\frac{1}{T} \sum_{i=1}^T \theta^{(t)}\). Questa approssimazione diventa sempre più accurata man mano che aumenta il numero di campioni $T, grazie alla Legge Forte dei Grandi Numeri.

Un altro vantaggio del Metodo di Monte Carlo è la sua capacità di approssimare la probabilità che una variabile casuale \(\theta\) cada all’interno di un intervallo specifico \((l, u)\). Questo può essere ottenuto calcolando la media campionaria della funzione indicatrice \(I(l < \theta < u)\) per ogni realizzazione \(\theta^{(t)}\), cioè \(Pr(l < \theta < u) \approx \frac{\text{numero di realizzazioni } \theta^{(t)} \in (l, u)}{T}\).

Nonostante la loro efficacia, un limite dei metodi Monte Carlo tradizionali risiede nella generazione efficiente di un elevato numero di campioni \(X_1, X_2, \ldots, X_n\). In risposta a questa sfida, i metodi di Monte Carlo basati su catene di Markov (MCMC) offrono una soluzione potente per simulare da distribuzioni complesse attraverso catene di Markov. L’evoluzione di questi algoritmi ha trasformato radicalmente la statistica e il calcolo scientifico, permettendo la simulazione da una vasta gamma di distribuzioni, anche in spazi di alta dimensionalità.

Le Catene di Markov

Le catene di Markov, ideate da Andrey Markov nel 1906, rappresentano un tentativo di estendere la legge dei grandi numeri a contesti in cui le variabili casuali non sono indipendenti. Tradizionalmente, la statistica si concentra su sequenze di variabili casuali indipendenti e identicamente distribuite (i.i.d.), simboleggiate come \(X_0, X_1, \ldots, X_n, \ldots\). In tali sequenze, ogni variabile è indipendente dalle altre e segue la stessa distribuzione, con \(n\) che rappresenta un indice temporale discreto. Tuttavia, questa assunzione di indipendenza non è sempre realistica nei modelli di fenomeni complessi, portando alla necessità di esplorare forme alternative di dipendenza tra variabili.

Per superare le limitazioni dell’indipendenza, le catene di Markov introducono una cosiddetta “dipendenza a un passo”, incarnata nella “proprietà di Markov”. Questa proprietà stabilisce che la previsione di un evento futuro \(X_{n+1}\) dipende unicamente dall’evento immediatamente precedente \(X_n\), indipendentemente dagli eventi passati \(X_0, X_1, X_2, \ldots, X_{n-1}\). La proprietà di Markov è espressa matematicamente come:

\[ P(X_{n+1} = j | X_n = i, X_{n-1} = i_{n-1}, \ldots, X_0 = i_0) = P(X_{n+1} = j | X_n = i). \]

Questa proprietà afferma che la previsione di un evento futuro dipende solo dall’evento immediatamente precedente, semplificando i calcoli relativi alle probabilità condizionali.

Le catene di Markov rappresentano un framework fondamentale per la modellazione delle dipendenze tra variabili casuali, una nozione cruciale in numerosi campi della statistica e della scienza dei dati, inclusa la metodologia MCMC (Markov Chain Monte Carlo). In particolare, nell’ambito dell’analisi bayesiana, l’uso di MCMC si rivela di estrema importanza, soprattutto quando non è possibile calcolare in modo analitico la distribuzione a posteriori.

L’algoritmo di Metropolis rappresenta una delle implementazioni più semplici del metodo MCMC. Questo algoritmo sfrutta la natura dipendente delle catene di Markov per navigare in modo efficace attraverso lo spazio della distribuzione a posteriori. Questo aspetto rende il MCMC uno strumento di grande potenza per affrontare problemi complessi in cui i metodi analitici tradizionali non possono essere applicati (per ulteriori dettagli, si veda Catene di Markov). In breve, il MCMC consente di generare un vasto insieme di valori per il parametro \(\theta\). Idealmente, questi valori riflettono la distribuzione a posteriori \(p(\theta \mid y)\) quando questa non può essere ottenuta direttamente. Questa caratteristica rende il MCMC uno strumento essenziale per risolvere problemi complessi in cui i metodi analitici convenzionali non sono applicabili.

Estrazione di campioni dalla distribuzione a posteriori

Nella discussione seguente ci porremo l’obiettivo di comprendere come utilizzare l’algoritmo di Metropolis per approssimare la distribuzione a posteriori \(p(\theta \mid y)\). A questo fine, il capitolo è strutturato in varie sezioni che facilitano la comprensione progressiva del tema.

• Inizieremo discutendo di come la distribuzione a posteriori possa essere approssimata mediante tecniche di simulazione convenzionali. Questa prima parte presuppone che la distribuzione target, o “a posteriori,” sia già conosciuta o disponibile per l’analisi.

• In seguito, passeremo a illustrare come l’algoritmo di Metropolis possa essere utilizzato per raggiungere lo stesso scopo—l’approssimazione della distribuzione a posteriori—mantenendo la presupposizione che la distribuzione target sia già nota.

• Concluderemo il capitolo esplorando le modalità con cui l’algoritmo di Metropolis può essere adattato per affrontare situazioni in cui la distribuzione a posteriori non è direttamente nota. In questi casi, spesso abbiamo a disposizione informazioni riguardanti la distribuzione a priori e la funzione di verosimiglianza, che possono essere utilizzate per ottenere un’approssimazione efficace della distribuzione a posteriori.

A titolo esemplificativo, utilizzeremo il dataset moma_sample.csv, il quale costituisce un campione casuale di 100 artisti provenienti dal Museo di Arte Moderna di New York (MoMA) e contiene diverse informazioni relative a ciascun artista.

Il nostro interesse è focalizzato sulla determinazione della probabilità che un artista presente nel MoMA appartenga alla generazione X o a una generazione successiva (nati dopo il 1965). Questa probabilità sarà indicata come \(\pi\). Iniziamo importando i dati.

moma_sample = pd.read_csv("../data/moma_sample.csv")

Esaminiamo le prime cinque righe del DataFrame.

moma_sample.head()

	artist	country	birth	death	alive	genx	gender	count	year_acquired_min	year_acquired_max
0	Ad Gerritsen	dutch	1940	2015.0	False	False	male	1	1981	1981
1	Kirstine Roepstorff	danish	1972	NaN	True	True	female	3	2005	2005
2	Lisa Baumgardner	american	1958	2015.0	False	False	female	2	2016	2016
3	David Bates	american	1952	NaN	True	False	male	1	2001	2001
4	Simon Levy	american	1946	NaN	True	False	male	1	2012	2012

Dai dati osserviamo che solo 14 artisti su 100 appartengono alla generazione X o a una generazione successiva.

result = moma_sample["genx"].value_counts()
print(result)

genx
False    86
True     14
Name: count, dtype: int64

Il valore campionato \(y = 14\) riflette le caratteristiche del campione che è stato osservato. Tuttavia, poiché il MOMA contiene opere di migliaia di artisti, sorge una domanda riguardante il vero valore di \(\theta\) (la probabilità di appartenere alla generazione X o a una generazione successiva) all’interno di questa popolazione.

Possiamo interpretare i dati \(y = 14\) come l’esito di una variabile casuale Binomiale con parametri \(N = 100\) e \(\theta\) sconosciuto.

Supponiamo che le nostre credenze pregresse riguardo a \(\theta\) possano essere modellate attraverso una distribuzione Beta(4, 6).

Sfruttando le proprietà delle distribuzioni coniugate, possiamo calcolare la distribuzione a posteriori esatta:

Y ~ Binomiale(100, π)
θ = Beta(4, 6)
θ | (Y = 14) ~ Beta(4 + 14, 6 + 100 - 14) → Beta(18, 92)

Nella figura successiva è rappresentata la distribuzione a posteriori del parametro \(\theta\), insieme alla distribuzione a priori scelta.

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

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

plt.fill_between(x, prior_density, alpha=0.5, label="Prior: Beta(4, 6)")
plt.fill_between(x, posterior_density, alpha=0.5, label="Posterior: Beta(18, 92)")
plt.xlabel("Parameter Value")
plt.ylabel("Density")
plt.legend()
plt.title("Prior and Posterior Densities")
plt.show()

Se vogliamo conoscere il valore della media a posteriori di \(\theta\), il risultato esatto è

\[ \bar{\theta}_{post} = \frac{\alpha}{\alpha + \beta} = \frac{18}{18 + 92} \approx 0.1636. \]

Simulazione con distribuzione target nota

Usiamo ora una simulazione numerica per stimare la media a posteriori di \(\theta\). Conoscendo la forma della distribuzione a posteriori \(Beta(18, 92)\), possiamo generare un campione di osservazioni casuali da questa distribuzione. Successivamente, calcoliamo la media delle osservazioni ottenute per ottenere un’approssimazione della media a posteriori.

Se vogliamo ottenere un risultato approssimato con poche osservazioni (ad esempio, 10), possiamo procedere con la seguente simulazione:

y = stats.beta(18, 92).rvs(10)
print(*y)

0.13247525044099537 0.17715535230368376 0.21962253612900964 0.24481835265576643 0.16679723516761935 0.11053999859705527 0.1163093395740192 0.14044864702085783 0.18260487197388198 0.12020234938139183

y.mean()

0.16109739332442805

Tuttavia, con solo 10 campioni l’approssimazione potrebbe non essere molto accurata. Più aumentiamo il numero di campioni (cioè il numero di osservazioni casuali generate), più precisa sarà l’approssimazione. Aumentando il numero di campioni, ad esempio a 10000, otteniamo un risultato più preciso:

stats.beta(18, 92).rvs(10000).mean()

0.16361874962751108

Quando il numero di campioni a posteriori diventa molto grande, la distribuzione campionaria converge alla densità della popolazione. Questo concetto si applica non solo alla media, ma anche ad altre statistiche descrittive, come la moda e la varianza.

È importante sottolineare che l’applicazione della simulazione di Monte Carlo è efficace per calcolare distribuzioni a posteriori solo quando conosciamo la distribuzione stessa e possiamo utilizzare funzioni Python per estrarre campioni casuali da tale distribuzione. Ciò è stato possibile nel caso della distribuzione a posteriori \(Beta(18, 92)\).

Tuttavia, questa situazione ideale non si verifica sempre nella pratica, poiché le distribuzioni a priori coniugate alla verosimiglianza sono spesso rare. Per esempio, nel caso di una verosimiglianza binomiale e una distribuzione a priori gaussiana, l’espressione

\[ p(\theta \mid y) = \frac{\mathrm{e}^{-(\theta - 1 / 2)^2} \theta^y (1 - \theta)^{n - y}} {\int_0^1 \mathrm{e}^{-(t - 1 / 2)^2} t^y (1 - t)^{n - y} dt} \]

rende impossibile calcolare analiticamente la distribuzione a posteriori di \(\theta\), precludendo quindi l’utilizzo diretto di Python per generare campioni casuali.

In queste circostanze, però, è possibile ottenere campioni casuali dalla distribuzione a posteriori mediante l’uso di metodi Monte Carlo basati su Catena di Markov (MCMC). Gli algoritmi MCMC, come ad esempio l’algoritmo Metropolis, costituiscono una classe di metodi che consentono di estrarre campioni casuali dalla distribuzione a posteriori senza richiedere la conoscenza della sua rappresentazione analitica. Le tecniche MCMC sono ampiamente adottate per risolvere problemi di inferenza bayesiana e rappresentano il principale strumento computazionale per ottenere stime approssimate di distribuzioni a posteriori in situazioni complesse e non analiticamente trattabili.

Algoritmo di Metropolis

L’algoritmo di Metropolis è un metodo avanzato per il campionamento da distribuzioni probabilistiche complesse. Appartiene alla famiglia dei metodi Monte Carlo Markov Chain (MCMC) e combina strategie di campionamento Monte Carlo con catene di Markov per navigare nello spazio dei parametri in modo intelligente. Questo consente di ottenere campioni rappresentativi della distribuzione di interesse indipendentemente dal punto di partenza, riducendo il rischio di bias nei risultati.

Passaggi Fondamentali dell’Algoritmo di Metropolis

1. Inizio: Si parte da un valore iniziale casuale per le variabili di interesse.

2. Generazione di Nuovi Punti: A ogni passo, si propone un nuovo valore (“salto”) basandosi su una distribuzione proposta, comunemente una distribuzione normale centrata sul punto corrente.

3. Valutazione: Si calcola la probabilità associata al nuovo punto utilizzando la funzione di densità di probabilità (pdf) della distribuzione che si intende campionare.

4. Decisione di Accettazione: Se il nuovo punto ha una probabilità maggiore rispetto al precedente (minore “energia”), lo si accetta. Questo favorisce la tendenza a spostarsi verso aree di maggiore probabilità.

5. Accettazione Probabilistica: Se invece il nuovo punto presenta una probabilità minore, l’accettazione avviene comunque con una certa probabilità, determinata dal rapporto tra le probabilità del nuovo e del vecchio punto. Questo consente di evitare locali ottimalità e di esplorare più ampiamente lo spazio dei parametri.

6. Raccolta dei Punti: I punti accettati vengono raccolti e costituiscono il campione dalla distribuzione desiderata.

La procedura continua con l’iterazione dei passi dal 2 al 5, permettendo alla catena di campioni di convergere alla distribuzione target. Inizialmente, i campioni potrebbero non essere rappresentativi, ma dopo un sufficiente numero di iterazioni (il cosiddetto “burn-in”), la distribuzione dei punti accettati dovrebbe avvicinarsi a quella che si intende esplorare.

L’efficacia di questo algoritmo risiede nella sua capacità di bilanciare esplorazione ed esploitazione attraverso l’accettazione probabilistica dei punti proposti, consentendo di ottenere un campionamento rappresentativo della distribuzione di probabilità complessa in esame.

Per una visualizzazione del comportamento dell’algoritmo di Metropolis nell’esplorare lo spazio dei parametri, si può consultare questo post. La distinzione tra i diversi metodi MCMC si basa principalmente sul tipo di distribuzione proposta e sul criterio di accettazione dei punti nel campionamento.

Versione di McElreath

McElreath [McE20] spiega l’algoritmo di Metropolis nel modo seguente. Re Markov era un autocrate benevolo di un regno insulare, un arcipelago circolare composto da 10 isole. Ogni isola era affiancata da altre due, e l’intero arcipelago formava un anello. Le isole erano di dimensioni diverse e quindi avevano popolazioni di diversa grandezza. La seconda isola aveva una popolazione circa il doppio di quella della prima, la terza circa il triplo della prima, e così via, fino all’isola più grande, che aveva una popolazione 10 volte maggiore della più piccola.

Il Buon Re aveva anche una serie di obblighi verso il suo popolo. Tra questi, Re Markov si impegnava a visitare ogni tanto ciascuna isola del suo regno. Poiché la gente amava il loro re, ogni isola preferiva che lui la visitasse più spesso. Così tutti concordarono che il re dovesse visitare ogni isola in proporzione alla sua dimensione demografica, visitando per esempio l’isola più grande 10 volte più spesso della più piccola.

Tuttavia, il Buon Re Markov non era amante delle agende o della contabilità, e quindi desiderava un modo per assolvere ai suoi obblighi senza pianificare i suoi viaggi con mesi di anticipo. Inoltre, poiché l’arcipelago era un anello, il Re insisteva nel muoversi solo tra isole adiacenti, per minimizzare il tempo trascorso sull’acqua.

Il consigliere del re, un certo Signor Metropolis, ideò una soluzione ingegnosa a queste esigenze. Chiameremo questa soluzione l’algoritmo di Metropolis.

1. Lancia una moneta per scegliere se andare all’isola a sinistra o a destra.

2. Determina la popolazione dell’isola proposta.

3. Determina la popolazione dell’isola corrente.

4. Spostati verso l’isola proposta con la probabilità calcolata.

Ripeti i passaggi dal 1 al 4.

A lungo termine, questo processo porterà a visitare ciascuna isola in proporzione alla sua popolazione. Se utilizziamo la popolazione dell’isola come metafora della densità di probabilità, possiamo impiegare lo stesso algoritmo per campionare da qualsiasi distribuzione di probabilità usando un generatore di numeri casuali.

Per l’inferenza bayesiana, applichiamo l’algoritmo di Metropolis assumendo:

• Le “isole” come i valori dei parametri.

• La “dimensione della popolazione” come le probabilità a posteriori.

In questa metafora, l’algoritmo di Metropolis ci aiuta a esplorare lo spazio dei parametri, muovendoci verso valori di parametri con maggior probabilità posteriore (simili a isole con maggior popolazione) e campionando da queste distribuzioni in maniera proporzionale alla loro probabilità.

L’algoritmo di Metropolis, nella versione di McElreath [McE20] è stato implementato da Dustin Stansbury come indicato di seguito.

def simulate_markov_visits(island_sizes, n_steps=100_000):
    """aka Metropolis algorithm"""
    
    # Metropolis algorithm
    island_idx = np.arange(len(island_sizes))
    visits = {ii: 0 for ii in island_idx}
    current_island_idx = np.random.choice(island_idx)
    markov_chain = []  # store history
    for _ in range(n_steps):
        
        # 1. Flip a coin to propose a direction, left or right
        coin_flip = np.random.rand() > 0.5
        direction = -1 if coin_flip else 1
        proposal_island_idx = np.roll(island_idx, direction)[current_island_idx]
        
        # 2. Proposal island size
        proposal_size = island_sizes[proposal_island_idx]
        
        # 3. Current island size
        current_size = island_sizes[current_island_idx]
        
        # 4. Go to proposal island if ratio of sizes is greater than another coin flip
        p_star = proposal_size / current_size
        move = np.random.rand() < p_star
        current_island_idx = proposal_island_idx if move else current_island_idx
        
        # 5. tally visits and repeat
        visits[current_island_idx] += 1
        markov_chain.append(current_island_idx)
        
    # Visualization
    island_idx = visits.keys()
    island_visit_density = [v / n_steps  for v in visits.values()]
    island_size_density = island_sizes / island_sizes.sum()
    
    _, axs = plt.subplots(1, 2, figsize=(12, 4))
    
    plt.sca(axs[0])
    plt.plot(island_sizes, lw=3, color='C0')
    plt.xlabel("Island Index")
    plt.ylabel("Island Population");
    plt.xticks(np.arange(len(island_sizes)));
    
    plt.sca(axs[1])
    plt.plot(island_idx, island_size_density, color='C0', lw=3, label='Population')
    plt.bar(island_idx, island_visit_density, color='k', width=.4, label='Visits')

    plt.xlabel("Island Index")
    plt.ylabel("Density");
    plt.xticks(np.arange(len(island_sizes)));

    plt.legend()
    return markov_chain

Verifichiamo l’algoritmo con una distribuzione approssimativamente gaussiana.

n_steps = 100_000
island_sizes = stats.norm(20, 10).pdf(np.linspace(1, 40, 10))
island_sizes /= island_sizes.min() / 1000

gaussian_markov_chain = simulate_markov_visits(island_sizes, n_steps) 

def plot_markov_chain(markov_chain, **hist_kwargs):
    _, axs = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
    plt.sca(axs[0])
    plt.plot(markov_chain[:500])
    plt.xlabel("visit number")
    plt.ylabel("island index")
    plt.title("Metropolis Algorithm Markov Chain\nFirst 500 Steps");
    
    plt.sca(axs[1])
    plt.hist(markov_chain, orientation='horizontal', density=True, **hist_kwargs);
    plt.title("Resulting Posterior\n(horizontal)");

plot_markov_chain(gaussian_markov_chain, bins=10)

Verifichiamo l’algoritmo con una distribuzione di Poisson.

island_sizes = stats.poisson(5).pmf(np.linspace(1, 15, 15))
island_sizes /= island_sizes.min() / 10

poisson_markov_chain = simulate_markov_visits(island_sizes, n_steps);

plot_markov_chain(poisson_markov_chain, bins=15)

Algoritmo di Metropolis con Distribuzione Target Nota

Torniamo al caso precedentemente esaminato in cui la distribuzione a posteriori è una Beta(18, 92). Ora, useremo una versione più generale dell’algoritmo di Metropolis rispetto a quello precedentemente discusso. Questa nuova versione ci consentirà di generare campioni da questa distribuzione nota in modo più flessibile.

Per iniziare, definiamo la distribuzione da cui verranno proposti i nuovi campioni. Nel nostro primo esempio, scegliamo una distribuzione proposta molto semplice: una piccola variazione, estratta da una distribuzione uniforme, verrà aggiunta al valore del campione precedente.

Ora procediamo a descrivere questa nuova implementazione dell’algoritmo di Metropolis.

# Proposal distribution
def uniprop(xprev, delta=0.05):
    return xprev + stats.uniform.rvs(loc=-delta, scale=2 * delta)

# Metropolis algorithm 
def metropolis_v1(p, qdraw, nsamp, xinit):
    """
    Implements the Metropolis algorithm for MCMC sampling.

    The function generates `nsamp` samples from a target distribution `p` using 
    the Metropolis algorithm. The algorithm starts with an initial state `xinit` 
    and iteratively applies a proposal mechanism governed by `qdraw` to explore 
    the target distribution.

    Args:
        p (callable): A function representing the target distribution to sample from.
                     Should accept a float or array-like and return the density 
                     (unnormalized is okay) at that point.
                     
        qdraw (callable): A function implementing the proposal mechanism. 
                          Takes the current state (float or array-like) as input 
                          and returns a proposed next state.
                          
        nsamp (int): The number of samples to draw from the target distribution.
        
        xinit (float or array-like): The initial state for the MCMC chain.

    Returns:
        numpy.ndarray: An array of shape (nsamp,) containing samples drawn from the 
        target distribution `p`.

    Example:
        >>> import numpy as np
        >>> from scipy import stats
        >>> target = stats.norm().pdf  # Standard Normal as target
        >>> qdraw = lambda x: x + np.random.normal(0, 1)  # Random-walk Metropolis
        >>> samples = metropolis_v1(target, qdraw, 10000, 0.0)
    """
    samples = np.empty(nsamp)
    x_prev = xinit
    for i in range(nsamp):
        x_star = qdraw(x_prev)
        p_star = p(x_star)
        p_prev = p(x_prev)
        pdfratio = p_star / p_prev
        if np.random.uniform() < min(1, pdfratio):
            samples[i] = x_star
            x_prev = x_star
        else:
            samples[i] = x_prev
    return samples

Sintassi della funzione

def metropolis_v1(p, qdraw, nsamp, xinit):

Parametri

1. p: Questa è la distribuzione target da cui si desidera campionare. È fornita come una funzione invocabile che prende un valore x e restituisce la densità di probabilità (o massa) in quel punto. La distribuzione target potrebbe anche non essere normalizzata.

2. qdraw: Questa è la distribuzione proposta da cui viene generato un possibile stato successivo (x_star), dato lo stato corrente (x_prev). Anche questa è una funzione invocabile, che dovrebbe restituire una nuova proposta basata sullo stato corrente.

3. nsamp: Questo è il numero di campioni che si desidera generare dalla distribuzione target. Determina quante iterazioni verranno eseguite dall’algoritmo.

4. xinit: Questo è lo stato iniziale della catena di Markov. L’algoritmo partirà da questo punto.

Valore di ritorno

• La funzione restituisce un array (samples) di dimensione nsamp, contenente i campioni estratti dalla distribuzione target.

Dettagli dell’algoritmo

1. Inizializzazione:

samples = np.empty(nsamp)
x_prev = xinit

• samples: Un array NumPy vuoto che sarà popolato con i campioni generati dall’algoritmo.

• x_prev: Lo stato corrente della catena di Markov. Inizializzato a xinit. Iniziamo dunque con un punto arbitrario nello spazio dei parametri. Questo primo valore della catena di Markov può essere scelto casualmente tra i valori possibili del parametro.

2. Il ciclo:

for i in range(nsamp):

• Si itera nsamp volte per generare nsamp campioni.

3. Fase di proposta:

x_star = qdraw(x_prev)

• Un nuovo stato (x_star) viene proposto sulla base dello stato corrente (x_prev) utilizzando la funzione di distribuzione proposta qdraw.

4. Calcolo della densità:

p_star = p(x_star)
p_prev = p(x_prev)

• Calcola le densità dei punti proposti (x_star) e correnti (x_prev) utilizzando la funzione di distribuzione target p.

5. Rapporto tra densità:

pdfratio = p_star / p_prev

• Viene calcolato il rapporto tra le densità (p_star / p_prev). Questo rapporto decide se accettare o meno lo stato proposto.

6. Fase di accettazione:

if np.random.uniform() < min(1, pdfratio):
    samples[i] = x_star
    x_prev = x_star
else:
    samples[i] = x_prev

• Viene estratto un numero casuale dalla distribuzione uniforme tra 0 e 1.

• Se questo numero è minore del minimo tra 1 e pdfratio, lo stato proposto viene accettato e x_prev viene aggiornato a x_star.

• Altrimenti, lo stato proposto viene respinto e lo stato corrente (x_prev) viene mantenuto per il ciclo successivo.

In sintesi, la decisione di accettare o rifiutare il parametro proposto x_star si basa sul valore di pdfratio. Se pdfratio è maggiore di 1, ciò indica che il parametro proposto è più probabile del parametro corrente x_prev; di conseguenza, x_star viene sempre accettato. Se, invece, pdfratio è minore di 1, il parametro proposto viene accettato con una probabilità che è proporzionale a pdfratio stesso. Per esempio, un pdfratio di 0.10 implica una probabilità del 10% di accettare x_star.

Questa strategia di accettazione ci consente di generare campioni che sono statisticamente rappresentativi della distribuzione target. Questo perché la probabilità di accettare un candidato è direttamente proporzionale alla sua densità della distribuzione che stiamo cercando di campionare.

Nella pratica, questo processo decisionale viene implementato confrontando pdfratio con un numero casuale \(u\), estratto da una distribuzione uniforme \(Unif(0, 1)\). Se pdfratio \(> u\), il parametro proposto x_star viene accettato, e la catena di Markov si sposta verso questo nuovo valore. In caso contrario, il valore corrente x_prev viene mantenuto come prossimo elemento nella catena.

Applichiamo ora l’algoritmo di Metropolis usando una Beta(18, 92) quale distribuzione target e una distribuzione uniforme quale distribuzione proposta.

# Initialization
init_state = 0.5
n_samples = 100000

# Using Beta(18, 92) as the target distribution
samps = metropolis_v1(lambda x: stats.beta.pdf(x, 18, 92), uniprop, n_samples, init_state)

Note

La ragione per cui è necessario scrivere lambda x: stats.beta.pdf(x, 18, 92) anziché stats.beta.pdf(x, 18, 92) è che la funzione metropolis_v1 si aspetta una funzione che può essere chiamata all’interno del codice di metropolis_v1. Scrivere direttamente stats.beta.pdf(x, 18, 92) come argomento causerà un errore di sintassi perché x non è definito in quel contesto. La funzione metropolis_v1 deve poter invocare p con differenti valori di x durante l’esecuzione. Quindi, ha bisogno di una “funzione” e non di un singolo valore. La funzione lambda lambda x: stats.beta.pdf(x, 18, 92) è invocabile: dato un x, restituirà il valore della PDF della distribuzione Beta in quel punto. Utilizzando la funzione lambda, fissiamo i parametri alpha = 18 e beta = 92 nella distribuzione Beta, rendendola una funzione di una sola variabile x, che è ciò che la funzione metropolis_v1 richiede.

In sostanza, la funzione lambda serve come un “incapsulamento” o “avvolgimento” per la funzione stats.beta.pdf, restringendola a specifici parametri e rendendola invocabile con una sola variabile, che è esattamente ciò che l’algoritmo di Metropolis necessita.

Esaminiamo i primi 20 valori ottenuti.

print(samps[0:20])

[0.5        0.48983531 0.48983531 0.49989737 0.47516176 0.44204386
 0.44204386 0.42783514 0.39304154 0.35385739 0.32065309 0.32065309
 0.29827246 0.29827246 0.29827246 0.29827246 0.26727104 0.27378934
 0.27378934 0.27378934]

Generiamo ora un istogramma di tutti i 100000 campioni prodotti dall’algoritmo di Metropolis, a cui sovrapponiamo la densità Beta(18, 92).

plt.hist(samps, bins=80, alpha=0.4, label="MCMC distribution", density=True)
x = np.linspace(0, 1, 200)
plt.plot(x, stats.beta.pdf(x, 18, 92), "C0", label="True distribution")
plt.legend()
plt.show()

Si noti che la procedura di Metropolis genera un insieme di campioni che approssima la distribuzione target desiderata. In altre parole, i campioni ottenuti seguono la distribuzione che stiamo cercando di esplorare, fornendo quindi un modo efficace per ottenere un campionamento casuale da questa distribuzione.

Ora esaminiamo una diversa funzione di proposta: una distribuzione gaussiana centrata attorno allo stato precedente, con una deviazione standard \(\sigma\). In questo caso, ogni nuovo punto candidato \(x^*\) sarà estratto da una gaussiana con media uguale allo stato corrente \(x_{\text{prev}}\) e deviazione standard \(\sigma\). Questo introduce un grado di esplorazione randomica attorno al punto attuale, permettendo all’algoritmo di navigare efficacemente nello spazio dei parametri. In pratica, questo significa che, se \(\sigma\) è piccola, il valore candidato \(x^*\) sarà simile al valore corrente \(x_{\text{prev}}\).

# Proposal distribution
def gaussprop(xprev, sigma=0.1):
    return xprev + stats.norm.rvs(scale=sigma)

Applichiamo nuovamente l’algorimo di Metropolis con questa nuova distribuzione proposta.

samps = metropolis_v1(lambda x: stats.beta.pdf(x, 18, 92), gaussprop, n_samples, init_state)
print(samps[0:20])

[0.5        0.5        0.5        0.50907441 0.43828701 0.43828701
 0.3692625  0.3692625  0.31822953 0.31822953 0.21435396 0.21435396
 0.21435396 0.21435396 0.18374932 0.18374932 0.18374932 0.22120697
 0.22120697 0.16270991]

Otteniamo anche in questo caso un campione casuale dalla distribuzione target.

plt.hist(samps, bins=80, alpha=0.4, label="MCMC distribution", density=True)
x = np.linspace(0, 1, 200)
plt.plot(x, stats.beta.pdf(x, 18, 92), "C0", label="True distribution")
plt.legend()
plt.show()

Algoritmo di Metropolis con distribuzione target incognita

Ora andiamo oltre: invece di presupporre che la distribuzione target sia già nota, cosa che si verifica quando utilizziamo distribuzioni coniugate, nel contesto bayesiano la distribuzione target che vogliamo campionare è in realtà il prodotto non normalizzato della verosimiglianza e della distribuzione a priori.

Iniziamo quindi col definire la distribuzione a priori. In questo esempio, la distribuzione a priori è una distribuzione Beta(2, 10) che è stata scelta solo per motivi didattici e non ha alcuna motivazione ulteriore.

def prior(p):
    alpha = 4
    beta = 6
    return stats.beta.pdf(p, alpha, beta)

Definiamo la verosimiglianza. Il problema presente richiede una verosimiglianza binomiale.

def likelihood(p):
    y = 14
    n = 100
    return stats.binom.pmf(y, n, p)

La distribuzione a posteriori non normalizzata è il prodotto della distribuzione a priori e della verosimiglianza. Si noti che, per il motivo spiegato prima, non è necessario normalizzare la distribuzione a posteriori.

def posterior(p):
    return likelihood(p) * prior(p)

Per la distribuzione proposta, faremo ricorso alla stessa funzione che abbiamo precedentemente definito. Tuttavia, apporteremo alcune modifiche all’algoritmo di Metropolis per adeguarlo a questo nuovo contesto.

Nell’implementazione di un algoritmo di Metropolis in ambito bayesiano, vi sono diversi aspetti da considerare attentamente. Ecco alcuni punti cruciali:

1. Simmetria della Distribuzione Proposta: È fondamentale che la distribuzione proposta sia simmetrica. Questa è una condizione necessaria per il funzionamento dell’algoritmo di Metropolis, ma non per quello di Metropolis-Hastings.

2. Valore Iniziale: Il valore iniziale scelto dovrebbe essere plausibile in base alla distribuzione a priori, in modo da facilitare la convergenza dell’algoritmo.

3. Probabilità Zero: Nel caso in cui la verosimiglianza o il prior assumano un valore di zero, il rapporto tra le densità di probabilità (pdfratio) all’interno dell’algoritmo di Metropolis risulterà indefinito. Pertanto, è importante garantire che il valore x_star, generato dalla distribuzione proposta, cada sempre entro i limiti del supporto sia del prior che della verosimiglianza.

Di seguito è presentato il codice, rivisto in base a queste considerazioni.

def metropolis_v2(p, qdraw, nsamp, xinit):
    samples = np.empty(nsamp)
    x_prev = xinit
    for i in range(nsamp):
        x_star = qdraw(x_prev)
        # Check that x_star is within the support of the prior and likelihood
        if 0 <= x_star <= 1:
            p_star = p(x_star)
            p_prev = p(x_prev)
            pdfratio = p_star / p_prev if p_prev > 0 else 1

            if np.random.uniform() < min(1, pdfratio):
                samples[i] = x_star
                x_prev = x_star
            else:
                samples[i] = x_prev
        else:
            samples[i] = x_prev  # reject automatically if outside support
    return samples

Le due versioni dell’algoritmo di Metropolis, metropolis_v1 e metropolis_v2, sono simili nella loro struttura di base, ma presentano alcune differenze chiave che riguardano la gestione dei casi speciali:

1. Controllo del supporto: Nella versione metropolis_v2, c’è un controllo esplicito che assicura che il valore proposto x_star cada all’interno del supporto della distribuzione priori e della verosimiglianza (if 0 <= x_star <= 1). Se x_star è fuori dal supporto, viene automaticamente rifiutato e il campione precedente x_prev viene conservato. Questo controllo è assente nella versione metropolis_v1.

2. Gestione delle probabilità zero: Nella versione metropolis_v2, viene effettuato un ulteriore controllo sulla probabilità precedente p_prev. Se p_prev è zero, pdfratio viene impostato automaticamente a 1 (pdfratio = p_star / p_prev if p_prev > 0 else 1). Questo gestisce i casi in cui potrebbe verificarsi una divisione per zero. Nella versione metropolis_v1, non vi è un controllo simile, e una divisione per zero potrebbe portare a comportamenti indesiderati.

In sintesi, metropolis_v2 è una versione più robusta dell’algoritmo che tiene conto sia del supporto della distribuzione target sia di eventuali valori zero nella densità di probabilità.

Si osservi un punto importante: nel calcolo di pdfratio, il rapporto tra la densità a posteriori del parametro proposto x_star e quella del parametro corrente x_prev, la costante di normalizzazione si annulla grazie alla regola di Bayes. Questo lascia nel rapporto solamente le componenti relative alla verosimiglianza e alla distribuzione a priori, che sono entrambe facilmente calcolabili. Matematicamente, questo si esprime come:

(62)\[ \begin{equation} \text{pdfratio} = \frac{p(\theta^* \mid y)}{p(\theta^{\text{prev}} \mid y)} = \frac{\frac{p(y \mid \theta^*) p(\theta^*)}{p(y)}}{\frac{p(y \mid \theta^{\text{prev}}) p(\theta^{\text{prev}})}{p(y)}} = \frac{p(y \mid \theta^*) p(\theta^*)}{p(y \mid \theta^{\text{prev}}) p(\theta^{\text{prev}})} \end{equation} \]

# Perform Metropolis sampling
samps = metropolis_v2(posterior, lambda x: gaussprop(x, sigma=0.05), n_samples, init_state)

Show code cell output Hide code cell output

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[29], line 2
      1 # Perform Metropolis sampling
----> 2 samps = metropolis_v2(posterior, lambda x: gaussprop(x, sigma=0.05), n_samples, init_state)

Cell In[28], line 8, in metropolis_v2(p, qdraw, nsamp, xinit)
      6 # Check that x_star is within the support of the prior and likelihood
      7 if 0 <= x_star <= 1:
----> 8     p_star = p(x_star)
      9     p_prev = p(x_prev)
     10     pdfratio = p_star / p_prev if p_prev > 0 else 1

Cell In[27], line 2, in posterior(p)
      1 def posterior(p):
----> 2     return likelihood(p) * prior(p)

Cell In[26], line 4, in likelihood(p)
      2 y = 14
      3 n = 100
----> 4 return stats.binom.pmf(y, n, p)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/scipy/stats/_distn_infrastructure.py:3377, in rv_discrete.pmf(self, k, *args, **kwds)
   3375 place(output, (1-cond0) + np.isnan(k), self.badvalue)
   3376 if np.any(cond):
-> 3377     goodargs = argsreduce(cond, *((k,)+args))
   3378     place(output, cond, np.clip(self._pmf(*goodargs), 0, 1))
   3379 if output.ndim == 0:

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/scipy/stats/_distn_infrastructure.py:601, in argsreduce(cond, *args)
    598 if not isinstance(newargs, list):
    599     newargs = [newargs, ]
--> 601 if np.all(cond):
    602     # broadcast arrays with cond
    603     *newargs, cond = np.broadcast_arrays(*newargs, cond)
    604     return [arg.ravel() for arg in newargs]

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/numpy/core/fromnumeric.py:2504, in all(a, axis, out, keepdims, where)
   2421 @array_function_dispatch(_all_dispatcher)
   2422 def all(a, axis=None, out=None, keepdims=np._NoValue, *, where=np._NoValue):
   2423     """
   2424     Test whether all array elements along a given axis evaluate to True.
   2425 
   (...)
   2502 
   2503     """
-> 2504     return _wrapreduction(a, np.logical_and, 'all', axis, None, out,
   2505                           keepdims=keepdims, where=where)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/numpy/core/fromnumeric.py:86, in _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs)
     84             return reduction(axis=axis, dtype=dtype, out=out, **passkwargs)
     85         else:
---> 86             return reduction(axis=axis, out=out, **passkwargs)
     88 return ufunc.reduce(obj, axis, dtype, out, **passkwargs)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/numpy/core/_methods.py:64, in _all(a, axis, dtype, out, keepdims, where)
     61 def _all(a, axis=None, dtype=None, out=None, keepdims=False, *, where=True):
     62     # Parsing keyword arguments is currently fairly slow, so avoid it for now
     63     if where is True:
---> 64         return umr_all(a, axis, dtype, out, keepdims)
     65     return umr_all(a, axis, dtype, out, keepdims, where=where)

KeyboardInterrupt: 

In somma, l’algoritmo Metropolis accetta come input il numero nsamp di passi da simulare, la deviazione standard \( \sigma \) della distribuzione proposta (in questo esempio, la deviazione standard \(\sigma\) è stata scelta empiricamente in modo tale da ottenere un tasso di accettazione adeguato), e la densità a priori. Come output, l’algoritmo restituisce una catena di valori del parametro, specificamente la sequenza \( \theta^{(1)}, \theta^{(2)}, \ldots, \theta^{\text{nsamp}} \). Uno degli aspetti cruciali per la riuscita dell’algoritmo è il raggiungimento della stazionarietà da parte della catena. In genere, i primi 1000-5000 valori vengono scartati in quanto rappresentano il periodo di “burn-in” della catena. Dopo un determinato numero di passi \( k \), la catena converge e i valori diventano campioni effettivi dalla distribuzione a posteriori \( p(\theta \mid y) \).

Il modello descritto è stato inizialmente proposto da Metropolis et al. nel 1953 Metropolis et al. [MRR+53]. Hastings nel 1970 introdusse un’estensione nota come algoritmo Metropolis-Hastings Hastings [Has70]. Altre varianti includono il campionatore di Gibbs, introdotto da Geman e Geman nel 1984 [GG84], l’Hamiltonian Monte Carlo Duane et al. [DKPR87], e il No-U-Turn Sampler (NUTS) utilizzato in pacchetti come PyMC e Stan Hoffman et al. [HG+14]. Per un’analisi più dettagliata e intuitiva dell’algoritmo Metropolis, si rimanda a Kruschke [Kru14].

Un elemento chiave da considerare nell’uso dell’algoritmo Metropolis è il tasso di accettazione, che è il rapporto tra il numero di valori del parametro proposti e il numero di quei valori che vengono effettivamente accettati. Un limite di questo algoritmo è la sua inefficienza relativa: rispetto alle sue varianti più moderne, l’algoritmo Metropolis tende ad essere meno efficiente.

Aspetti computazionali

Warm-up/Burn-in

Una catena di Markov ha bisogno di alcune iterazioni per raggiungere la distribuzione stazionaria. Queste iterazioni sono comunemente chiamate iterazioni di warm-up o burn-in (a seconda dell’algoritmo e del software utilizzato) e vengono di solito scartate. In molti programmi software, la prima metà delle iterazioni viene considerata come iterazioni di warm-up, quindi, nel caso presente, anche se abbiamo ottenuto 100000 iterazioni, ne utilizzeremo solo 50000.

Sintesi della distribuzione a posteriori

L’array samps contiene 100000 valori di una catena di Markov. Escludiamo i primi 50000 valori considerati come burn-in e consideriamo i restanti 50000 valori come un campione casuale estratto dalla distribuzione a posteriori \(p(\theta \mid y)\).

Mediante i valori della catena così ottenuta è facile trovare una stima a posteriori del parametro \(\theta\). Per esempio, possiamo trovare la stima della media a posteriori.

burnin = int(n_samples * 0.5)
burnin

50000

np.mean(samps[burnin:])

0.16314416036984744

Oppure possiamo stimare la deviazione standard della distribuzione a posteriori.

np.std(samps[burnin:])

0.03527332289207083

Visualizziamo un trace plot dei valori della catena di Markov dopo il periodo di burn-in.

plt.plot(samps[burnin:])
plt.xlabel("sample")
plt.ylabel("theta")
plt.show()

Il trace plot indica che la catena inizia con un valore casuale per poi spostarsi rapidamente nell’area intorno a 0.16, che è l’area con alta densità a posteriori. Successivamente, oscilla intorno a quel valore per le iterazioni successive.

plt.plot(samps[:500])
plt.xlabel("sample")
plt.ylabel("theta")
plt.show()

L’istogramma mostrato di seguito, sul quale è stata sovrapposta la distribuzione a posteriori derivata analiticamente – specificamente una \(\text{Beta}(25, 17)\) – dimostra che la catena converge effettivamente alla distribuzione a posteriori desiderata.

plt.hist(samps[burnin:], bins=80, alpha=0.4, label="MCMC distribution", density=True)
# plot the true function
x = np.linspace(0, 1, 100)
plt.plot(x, stats.beta.pdf(x, 18, 92), "C0", label="True distribution")
plt.legend()
plt.show()

È possibile usare la funzione summary del pacchetto AriviZ per calolare l’intervallo di credibilità, ovvero l’intervallo che contiene la proporzione indicata dei valori estratti a caso dalla distribuzione a posteriori.

az.summary(samps[burnin:], kind="stats", hdi_prob=0.94, round_to=2)

	mean	sd	hdi_3%	hdi_97%
x	0.16	0.04	0.1	0.23

Un KDE plot corrispondente all’istogramma precedente si può generare usando az.plot_posterior(). La curva rappresenta l’intera distribuzione a posteriori e viene calcolata utilizzando la stima della densità del kernel (KDE) che serve a “lisciare” l’istogramma.

az.plot_posterior(samps[burnin:])
plt.show()

L’HDI è una scelta comune nelle statistiche bayesiane e valori arrotondati come 50% o 95% sono molto popolari. Ma ArviZ utilizza il 94% come valore predefinito. La ragione di questa scelta è che il 94% è vicino al valore ampiamente utilizzato del 95%, ma è anche diverso da questo, così da servire da “amichevole promemoria” che non c’è niente di speciale nella soglia del 5%. Idealmente sarebbe opportuno scegliere un valore che si adatti alle specifiche esigenze dell’analisi statistica che si sta svolgendo, o almeno riconoscere che si sta usando un valore arbitrario.

Diagnostiche della soluzione MCMC

Catene multiple

Poiché ciascuno stato in una catena di Markov dipende dagli stati precedenti, il valore o i valori iniziali possono influenzare i valori campionati. Una soluzione per verificare la sensibilità rispetto ai valori iniziali è utilizzare più catene, ognuna con diversi valori iniziali. Se più catene campionano la stessa distribuzione target, queste dovrebbero mescolarsi bene, ovvero intersecarsi l’una con l’altra in un trace plot.

Stazionarietà

Un punto importante da verificare è se il campionatore ha raggiunto la sua distribuzione stazionaria. La convergenza di una catena di Markov alla distribuzione stazionaria viene detta “mixing”.

Autocorrelazione

Ogni passo nell’algoritmo MCMC è chiamato iterazione. I valori campionati sono dipendenti, il che significa che il valore all’iterazione \(m\) dipende dal valore all’iterazione \(m-1\). Questa è una differenza importante rispetto alle funzioni che generano campioni casuali indipendenti, come beta(25, 17).rvs(). I valori campionati formano una catena di Markov, il che significa che ciascun valore campionato è correlato con il valore precedente (ad esempio, se \(\theta(m)\) è grande, \(\theta(m+1)\) sarà anch’esso grande).

Informazioni sul “mixing” della catena di Markov sono fornite dall’autocorrelazione. L’autocorrelazione misura la correlazione tra i valori successivi di una catena di Markov. Il valore \(m\)-esimo della serie ordinata viene confrontato con un altro valore ritardato di una quantità \(k\) (dove \(k\) è l’entità del ritardo) per verificare quanto si correli al variare di \(k\). L’autocorrelazione di ordine 1 (lag 1) misura la correlazione tra valori successivi della catena di Markow (cioè, la correlazione tra \(\theta^{(m)}\) e \(\theta^{(m-1)}\)); l’autocorrelazione di ordine 2 (lag 2) misura la correlazione tra valori della catena di Markow separati da due “passi” (cioè, la correlazione tra \(\theta^{(m)}\) e \(\theta^{(m-2)}\)); e così via.

L’autocorrelazione di ordine \(k\) è data da \(\rho_k\) e può essere stimata come:

(63)\[\begin{split} \begin{align} \rho_k &= \frac{Cov(\theta_m, \theta_{m+k})}{Var(\theta_m)}\notag\\ &= \frac{\sum_{m=1}^{n-k}(\theta_m - \bar{\theta})(\theta_{m-k} - \bar{\theta})}{\sum_{m=1}^{n-k}(\theta_m - \bar{\theta})^2} \qquad\text{con }\quad \bar{\theta} = \frac{1}{n}\sum_{m=1}^{n}\theta_m. \end{align} \end{split}\]

Per fare un esempio pratico, simuliamo dei dati autocorrelati.

x = pd.array([22, 24, 25, 25, 28, 29, 34, 37, 40, 44, 51, 48, 47, 50, 51])
print(*x)

22 24 25 25 28 29 34 37 40 44 51 48 47 50 51

L’autocorrelazione di ordine 1 è semplicemente la correlazione tra ciascun elemento e quello successivo nella sequenza.

sm.tsa.acf(x)

array([ 1.        ,  0.83174224,  0.65632458,  0.49105012,  0.27863962,
        0.03102625, -0.16527446, -0.30369928, -0.40095465, -0.45823389,
       -0.45047733, -0.36933174])

Nell’esempio, il vettore x è una sequenza temporale di 15 elementi. Il vettore \(x'\) include gli elementi con gli indici da 0 a 13 nella sequenza originaria, mentre il vettore \(x''\) include gli elementi 1:14. Gli elementi delle coppie ordinate dei due vettori avranno dunque gli indici \((0, 1)\), \((1, 2), (2, 3), \dots (13, 14)\) degli elementi della sequenza originaria. La correlazione di Pearson tra i vettori \(x'\) e \(x''\) corrisponde all’autocorrelazione di ordine 1 della serie temporale.

Nell’output precedente

• 0.83174224 è l’autocorrelazione di ordine 1 (lag = 1),

• 0.65632458 è l’autocorrelazione di ordine 2 (lag = 2),

• 0.49105012 è l’autocorrelazione di ordine 3 (lag = 3),

• ecc.

È possibile specificare il numero di ritardi (lag) da utilizzare con l’argomento nlags:

sm.tsa.acf(x, nlags=4)

array([1.        , 0.83174224, 0.65632458, 0.49105012, 0.27863962])

In Python possiamo creare un grafico della funzione di autocorrelazione (correlogramma) per una serie temporale usando la funzione tsaplots.plot_acf() dalla libreria statsmodels.

tsaplots.plot_acf(x, lags=9)
plt.show()

Per i dati dell’esempio in discussione otteniamo la situazione seguente.

tsaplots.plot_acf(samps[burnin:], lags=9)
plt.show()

Il correlogramma è uno strumento grafico usato per la valutazione della tendenza di una catena di Markov nel tempo. Il correlogramma si costruisce a partire dall’autocorrelazione \(\rho_k\) di una catena di Markov in funzione del ritardo \(k\) con cui l’autocorrelazione è calcolata: nel grafico ogni barretta verticale riporta il valore dell’autocorrelazione (sull’asse delle ordinate) in funzione del ritardo (sull’asse delle ascisse).

In situazioni ottimali l’autocorrelazione diminuisce rapidamente ed è effettivamente pari a 0 per piccoli lag. Ciò indica che i valori della catena di Markov che si trovano a più di soli pochi passi di distanza gli uni dagli altri non risultano associati tra loro, il che fornisce una conferma del “mixing” della catena di Markov, ossia della convergenza alla distribuzione stazionaria. Nelle analisi bayesiane, una delle strategie che consentono di ridurre l’autocorrelazione è quella di assottigliare l’output immagazzinando solo ogni \(m\)-esimo punto dopo il periodo di burn-in. Una tale strategia va sotto il nome di thinning.

Nel seguente correlogramma, analizziamo la medesima catena di Markov. Tuttavia, in questa occasione applichiamo un “thinning” (sottocampionamento) con un fattore di 5.

thin=5
sampsthin=samps[burnin::thin]
tsaplots.plot_acf(sampsthin, lags=9)
plt.show()

Si può notare come l’autocorrelazione diminuisce molto più rapidamente.

Tasso di accettazione

Quando si utilizza l’algoritmo Metropolis, è importante monitorare il tasso di accettazione e assicurarsi che sia nell’intervallo ottimale. Se si accetta quasi sempre il candidato proposto, probabilmente significa che, in ogni iterazione, la catena salta solo di un piccolo passo (in modo che il rapporto di accettazione sia vicino a 1 ogni volta). Di conseguenza, la catena impiegherà molte iterazioni per raggiungere altre regioni della distribuzione stazionaria e i campioni consecutivi saranno molto fortemente correlati. D’altra parte, se il tasso di accettazione è molto basso, la catena rimarrà bloccata nella stessa posizione per molte iterazioni prima di spostarsi verso uno stato diverso. Per l’algoritmo Metropolis base con un singolo parametro con una distribuzione proposta Gaussiana normale, un tasso di accettazione ottimale è compreso tra il 40% e il 50%.

Test di convergenza

Un test di convergenza può essere svolto in maniera grafica mediante le tracce delle serie temporali (trace plot), cioè il grafico dei valori simulati rispetto al numero di iterazioni. Se la catena è in uno stato stazionario le tracce mostrano assenza di periodicità nel tempo e ampiezza costante, senza tendenze visibili o andamenti degni di nota.

Ci sono inoltre alcuni test che permettono di verificare la stazionarietà del campionatore dopo un dato punto. Uno è il test di Geweke che suddivide il campione, dopo aver rimosso un periodo di burn in, in due parti. Se la catena è in uno stato stazionario, le medie dei due campioni dovrebbero essere uguali. Un test modificato, chiamato Geweke z-score, utilizza un test \(z\) per confrontare i due subcampioni ed il risultante test statistico, se ad esempio è più alto di 2, indica che la media della serie sta ancora muovendosi da un punto ad un altro e quindi è necessario un periodo di burn-in più lungo.

Effective sample size (ESS)

Quando le iterazioni sono dipendenti, ogni iterazione contiene informazioni sovrapposte con le iterazioni precedenti. In altre parole, quando si ottengono 500 campioni dipendenti dalla distribuzione a posteriori, questi contengono solo informazioni equivalenti a < 500 campioni indipendenti. L’ESS (Effective Sample Size) quantifica la quantità effettiva di informazioni, quindi una catena con ESS = n conterrà approssimativamente la stessa quantità di informazioni di n campioni indipendenti. In generale, vogliamo che l’ESS sia almeno 400 per un’utilizzazione generale nel riassumere la distribuzione a posteriori.

Commenti e considerazioni finali

In generale, la distribuzione a posteriori dei parametri di un modello statistico non può essere determinata per via analitica. Tale problema viene invece affrontato facendo ricorso ad una classe di algoritmi per il campionamento da distribuzioni di probabilità che sono estremamente onerosi dal punto di vista computazionale e che possono essere utilizzati nelle applicazioni pratiche solo grazie alla potenza di calcolo dei moderni computer. Lo sviluppo di software che rendono sempre più semplice l’uso dei metodi MCMC, insieme all’incremento della potenza di calcolo dei computer, ha contribuito a rendere sempre più popolare il metodo dell’inferenza bayesiana che, in questo modo, può essere estesa a problemi di qualunque grado di complessità.

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

The watermark extension is already loaded. To reload it, use:
  %reload_ext watermark
Last updated: Fri Jan 26 2024

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

scipy      : 1.11.4
matplotlib : 3.8.2
numpy      : 1.26.2
pandas     : 2.1.4
arviz      : 0.17.0
statsmodels: 0.14.1
seaborn    : 0.13.0

Watermark: 2.4.3