Campionamento di Gibbs

Da Wikipedia, l'enciclopedia libera.

In statistica e in fisica statistica, un campionamento di Gibbs o un campionatore di Gibbs è un algoritmo di catena di Markov Monte Carlo (MCMC) per ottenere una sequenza di campioni casuali da una distribuzione di probabilità multivariata (cioè dalla distribuzione di probabilità congiunta di due o più variabili casuali) quando il campionamento diretto si dimostra difficoltoso. Questa sequenza può essere usata per approssimare la distribuzione congiunta (ad esempio per generare un istogramma della distribuzione); per approssimare la distribuzione marginale di una delle variabili, o di vari sottoinsiemi delle variabili (per esempio, parametri sconosciuti oppure variabili latenti); oppura ancora per calcolare un integrale (come il valore atteso di una delle variabili).

Il campionamento di Gibbs è comunemente usato come uno strumento per eseguire inferenza statistica, specialmente inferenza bayesiana. Per sua natura è un algoritmo casuale (cioè un algoritmo che fa uso di numeri casuali, e quindi può produrre risultati distinti ogni volta che viene eseguito), ad è un'alternativa agli algoritmi deterministici impiegati nell'inferenza statistica come l'algoritmo variazionale di Bayes o l'algoritmo di massimizzazione del valore atteso (expectation-maximization algorithm (EM)).

Similmente ad altri algoritmi MCMC, il campionamento di Gibbs genera catene di Markov di campioni, ognuno dei quali è autocorrelato a quelli generati immediatamente prima e dopo di lui. Conseguentemente si deve prestare attenzione se si vogliono campioni indipendenti (ad esempio considerando come buoni solo quelli generati ogni un certo numero, ad esempio ogni cento). Inoltre (di nuovo, come in altri algoritmi MCMC), campioni provenienti dalla parte iniziale della catena (il periodo cosiddetto di burn-in) possono non rappresentare accuratamente la distribuzione desiderata e vanno perciò scartati.

Introduzione[modifica | modifica sorgente]

Il campionamento di Gibbs prende nome dal fisico J. W. Gibbs, in relazione all'analogia tra l'algoritmo di campionamento e la fisica statistica. L'algoritmo fu descritto dai fratelli Stuart e Donald Geman nel 1984, più di otto decadi dalla morte di Gibbs.[1]

Nella sua versione base il campionamento di Gibbs è un caso speciale dell'algoritmo di Metropolis-Hastings. Tuttavia, nelle sue versioni estese (cfr. più sotto), può essere inglobato in un contesto di riferimento generale per il campionamento di grandi insiemi di variabili mediante il campionamento di singole variabili (o in vari casi, di gruppi di variabili) a turno, e può incorporare l'algoritmo di Metropolis-Hastings (o metodi simili come il campionamento a strati) per implementare uno o più passi di campionamento.

Il campionamento di Gibbs è applicabile quando la distribuzione congiunta non è nota esplicitamente oppure è difficile da campionare direttamente, ma è nota la distribuzione condizionata di ogni variabile e che si presuppone essere più facile da campionare. Il campionamento di Gibbs genera a turno per ogni variabile un campione dalla loro distribuzione condizionata sui valori correnti delle restanti variabili. Può essere dimostrato (cfr. Gelman et al., 1995) che la sequenza di campioni costituisce una catena di Markov e che la distribuzione stazionaria di tale catena è di fatto la distribuzione congiunta cercata.

Il campionamento di Gibbs è particolarmente adatto per campionare la distribuzione a posteriori di una rete bayesiana in quanto tale genere di rete è specificata tipicamente mediante una collezione di distribuzioni condizionate.

Implementazione[modifica | modifica sorgente]

Il campionamento di Gibbs, nella sua incarnazione più semplice, è un caso speciale dell'algoritmo di Metropolis-Hastings. Il punto chiave del campionamento di Gibbs è che data una distribuzione multivariata, è più semplice campionare da una distribuzione condizionata che da una marginale mediante integrazione sopra una distribuzione congiunta. Supponiamo di voler ottenere \left.k\right. campioni di \mathbf{X} = \{x_1, \dots, x_n\} da una distribuzione congiunta \left.p(x_1, \dots, x_n)\right.. Denotiamo l'elemento i-esimo come \mathbf{X}^{(i)} = \{x_1^{(i)}, \dots, x_n^{(i)}\} e procediamo come segue:

  1. Cominciamo col fissare il valore \mathbf{X}^{(0)} per ogni variabile.
  2. Per ogni campione i = \{1 \dots k\}, campioniamo ogni variabile x_j^{(i)} dalla distribuzione condizionata p(x_j^{(i)}|x_1^{(i)},\dots,x_{j-1}^{(i)},x_{j+1}^{(i-1)},\dots,x_n^{(i-1)}). Campioniamo cioè ogni variabile della distribuzione da quest'ultima ma condizionandola su tutte le altre variabili, facendo uso dei loro valori più recenti ed aggiornando la variabile al suo nuovo valore non appena questo è stato campionato.

I campioni in tal modo approssimano la distribuzione congiunta di tutte le variabili. Inoltre, la distribuzione marginale di qualsiasi sottoinsieme di variabili può essere approssimata semplicemente esaminando i campioni di quel sottoinsieme di variabili, ignorando le restanti. In aggiunta, il valore atteso di qualsiasi variabile può essere approssimato mediando sopra tutti i campioni.

  • I valori iniziali delle variabili possono essere determinati casualmente o mediante qualche altro algoritmo come quello di expectation-maximization.
  • Non è di fatto necessario determinare un valore iniziale per la prima variabile campionata.
  • È pratica comune ignorare i valori campionati nella parte iniziale (il cosiddetto periodo di burn-in ossia di "riscaldamento") e in seguito prendere in considerazione solo un campione ogni n quando vengono richieste delle medie dei valori attesi. Per esempio, i primi 1000 campioni possono essere ignorati e successivamente utilizzarne solo un valore ogni 100 per il calcolo ad esempio di una media, scartando tutti gli altri. I motivi per procedere in questo modo sono che (1) i campioni consecutivi non sono indipendenti tra loro ma formano una catena di Markov con un certo grado di correlazione; (2) la distribuzione stazionaria della catena di Markov è la distribuzione congiunta sopra le variabili desiderata, ma può richiedere tempo affinché venga raggiunta la stazionarietà. Talvolta, per determinare il grado di autocorrelazione tra i campioni e il valore di n (il valore del periodo dei campioni conservati per i calcoli) sono richiesti degli algoritmi aggiuntivi e comunque è necessaria una certa pratica.
  • Il procedimento di Simulated Annealing viene spesso usato per ridurre l'effetto di random walk nella parte iniziale del processo di campionamento (cioè la tendenza a vagare lentamente nello spazio di campionamento, con un elevato grado di autocorrelazione tra i campioni, piuttosto che muoversi rapidamente, come desiderato). Altre tecniche che possono ridurre l'autocorrelazione sono il campionamento di Gibbs collassato, il campionamento di Gibbs bloccato e il sovrarilassamento ordinato; cfr. più sotto.

La relazione tra distribuzione condizionata e distribuzione congiunta[modifica | modifica sorgente]

In aggiunta a quando detto finora, la distribuzione condizionata di una variabile date tutte le altre è proporzionale alla distribuzione congiunta:

p(x_j|x_1,\dots,x_{j-1},x_{j+1},\dots,x_n) = \frac{p(x_1,\dots,x_n)}{p(x_1,\dots,x_{j-1},x_{j+1},\dots,x_n)} \propto p(x_1,\dots,x_n)

Il segno di proporzionalità in questo caso significa che il denominatore non è una funzione di x_j e perciò è il medesimo per tutti i valori di x_j. In pratica, per determinare la natura della distribuzione condizionata del fattore x_j è più semplice fattorizzare la distribuzione congiunta in accordo alle distribuzioni condizionate individuali definite tramite il modello grafico delle variabili (ossia il grafo associato alla struttura di indipendenza condizionata tra le varie variabili casuali in gioco, cfr. anche [grafo aleatorio]), ignorare tutti i fattori che non sono funzioni di x_j (i quali tutti assieme, ed assieme con il denominatore sopra definito, costituiscono la costante di normalizzazione), e quindi alla fine reimpostare la costante di normalizzazione secondo la necessità. In pratica, questo significa fare una delle tre scelte:

  1. Se la distribuzione è discreta, vengono calcolate le probabilità individuali di tutti i possibili valori di x_j, e quindi sommate per trovare la costante di normalizzazione.
  2. Se la distribuzione è continua e di una forma nota, la costante di normalizzazione è di conseguenza nota anch'essa.
  3. Negli altri casi, siccome molti metodi di campionamento non la richiedono, solitamente la costante di normalizzazione può essere ignorata.

Inferenza[modifica | modifica sorgente]

Il campionamento di Gibbs è comunemente usato per inferenza statistica (ad esempio per determinare il migliore valore di un parametro, come il numero di persone che probabilmente faranno acquisti presso un certo magazzino in un dato giorno, il candidato preferito, ecc.). L'idea è di incorporare i dati osservati nel processo di campionamento mediante la creazione di variabili separate per ogni parte dei dati osservati fissando le variabili in questione ai loro valori osservati, piuttosto che campionare da quelle variabili. La distribuzione delle variabili rimanenti è allora a tutti gli effetti una distribuzione a posteriori condizionata sui dati osservati.

Il valore con probabilità maggiore di un parametro desiderato (la moda) potrebbe allora essere selezionato semplicemente come il valore del campione che occorre più frequentemente; questo è essenzialmente equivalente a massimizzare a posteriori la stima di un parametro. (Poiché i parametri sono solitamente continui, è spesso necessario raggruppare i valori campionati che sono adiacenti tra loro in un numero finito di intervalli consecutivi allo scopo di ottenere una stima sensata della moda). Più comunemente, tuttavia, è scelto il valore atteso (media) dei valori campionati; questo è uno stimatore bayesiano che ha il vantaggio di prendere in considerazione l'intera distribuzione disponibile dal campionamento bayesiano, mentre un algoritmo di massimizzazione come quello di massimizzazione del valore atteso è in grado di fornire solo un singolo punto dalla distribuzione. Per esempio, per una distribuzione unimodale il valore della media (valore atteso) è solitamente simile a quello della moda (il valore più comune), ma se la distribuzione è asimmetrica in una direzione, la media si sposterà correttamente in quella direzione in modo da tenere conto della maggiore probabilità in quella direzione. (Si noti tuttavia che se la distribuzione è multimodale, il valore atteso può coincidere con un valore avente probabilità nulla rendendo la sua interpretazione più ardua rispetto ad esempio a quella della moda nel medesimo contesto).

Nonostante alcune delle variabili corrispondano tipicamente a parametri di interesse, altre sono variabili non interessanti ("prive di senso") vengono introdotte in generale nel modello allo scopo di esprimere propriamente le relazioni tra variabili. Nonostante i valori campionati rappresentino la distribuzione congiunta sopra tutte le variabili, i parametri non interessanti possono essere semplicemente ignorati nel momento in cui vengono calcolati il valore atteso o la moda delle variabili interessanti; questo equivale a marginalizzare sopra le variabili non interessanti. Quando un valore è desiderato per più variabili, il valore atteso è calcolato semplicemente sopra ogni variabile presa separatamente. (Tuttavia, quando si calcola la moda tutte le variabili devono essere considerate).

L'apprendimento supervisionato, l'apprendimento non supervisionato e l'apprendimento parzialmente supervisionato (vale a dire l'apprendimento con valori mancanti) possono tutti essere trattati semplicemente fissando tutte le variabili i cui valori sono noti e campionando dalle rimanenti.

Per i dati osservati, ci sarà una variabile per ogni osservazione piuttosto che, per esempio, una variabile corrispondente alla media campionaria o alla varianza campionaria di un insieme di osservazioni. Infatti, non ci saranno di fatto variabili corrispondenti a concetti come "media campionaria" o "varianza campionaria". Invece, in tale caso ci saranno delle variabili rappresentanti i valori incogniti della media e della varianza, e la determinazione dei valori campione per quelle variabili risulterà automaticamente come risultato dell'applicazione del campionamento di Gibbs.

Modelli lineari generalizzati (cioè varianti della regressione lineare) possono di fatto essere trattati mediante il campionamento di Gibbs. Per esempio, la regressione probit utile per determinare la probabilità di una data scelta binaria (ad esempio si/no), con distribuzioni a priori di tipo normale dei coefficienti di regressione, può essere implementata con campionamenti di Gibbs in quanto è possibile introdurre variabili aggiuntive e sfruttare la coniugazione con la distribuzione a priori. Tuttavia, la regressione logistica non può essere trattata in questo modo. Una possibilità è quella di approssimare le funzione logistica con una miscela di (tipicamente 7-9) distribuzioni normali. Più comunemente, tuttavia, al posto del campionamento di Gibbs viene usato l'algoritmo di Metropolis-Hastings.

Fondamento matematico[modifica | modifica sorgente]

Supponiamo che un campione \left.X\right. sia preso da una distribuzione dipendente da un vettore di parametri \theta \in \Theta \,\! di lunghezza \left.d\right., con distribuzione a priori g(\theta_1, \ldots , \theta_d). Può accadere che \left.d\right. sia veramente esteso e che l'integrazione numerica per trovare le densità marginali \left.\theta_i\right. sia computazionalmente dispendiosa. Allora un metodo alternativo per calcolare le densità marginali è quello di creare una catena di Markov sullo spazio \left.\Theta_i\right. ripetendo i due passi seguenti:

  1. Scegli un indice casuale 1 \leq j \leq d
  2. Scegli un nuovo valore \left.\theta_j\right. conformemente a g(\theta_1, \ldots , \theta_{j-1} , \, \cdot \, , \theta_{j+1} , \ldots , \theta_d )

Questi passi definiscono una catena di Markov reversibile con la distribuzione invariante desiderata \left.g\right.. Questo può essere provato nel modo seguente. Definiamo x \sim_j y se \left.x_i = y_i\right. per tutti gli i \neq j e sia \left.p_{xy}\right. la probabilità di un salto da x \in \Theta a y \in \Theta. Allora, le probabilità di transizione sono

p_{xy} = \begin{cases}
\frac{1}{d}\frac{g(y)}{\sum_{z \in \Theta: z \sim_j x} g(z) } & x \sim_j y \\
0 & \text{altrimenti}
\end{cases}

In tal modo


g(x) p_{xy} = \frac{1}{d}\frac{ g(x) g(y)}{\sum_{z \in \Theta: z \sim_j x} g(z) }
= \frac{1}{d}\frac{ g(y) g(x)}{\sum_{z \in \Theta: z \sim_j y} g(z) }
= g(y) p_{yx}

in quanto x \sim_j y è una relazione di equivalenza. Perciò le equazioni di bilancio dettagliato sono soddisfatte, implicando che la catena è reversibile e che possiede una distribuzione invariante \left.g\right..

In pratica, il suffisso \left.j\right. non è scelto a caso e la catena è ciclica nei suffissi posti in un certo ordine. In generale questo equivale ad un processo di Markov non stazionario, ma ogni singolo passo sarà ancora reversibile e il processo complessivamente avrà ancora la distribuzione stazionaria desiderata (fintanto che la catena può accedere a tutti gli stati nell'ordine prefissato).

Varianti ed estensioni[modifica | modifica sorgente]

Esistono numerose varianti del campionamento base di Gibbs. L'obiettivo di queste varianti è ridurre l'autocorrelazione tra i campioni in modo sufficiente da giustificare il costo associato all'elaborazione aggiuntiva.

Campionamento di Gibbs bloccato[modifica | modifica sorgente]

Campionamento di Gibbs collassato[modifica | modifica sorgente]

  • Un campionamento di Gibbs collassato integra marginalizzando sopra una o più variabili quando esegue il campionamento per altre variabili. Per esempio, immaginiamo che un modello consista di tre variabili A, B e C. Un semplice campionamento di Gibbs campionerebbe p(A|B,C), da p(B|A,C) e da p(C|A,B). Un campionamento di Gibbs collassato può sostituire il passo di campionamento per A con un campione preso dalla distribuzione marginale p(A|C), con in questo caso la variabile B integrata al di fuori (dell'ambito di integrazione delle altre due variabili). In alternativa, la variabile B può essere interamente fuori collassata ovvero esclusa dalla fase di calcolo, campionando da p(A|C) e da p(C|A) e non campionando del tutto da B. La distribuzione di una variabile A che è generata quando di collassa una variabile B che parametrizza la distribuzione di A è chiamata distribuzione composta; il campionamento da questa distribuzione è generalmente trattabile quando B è la variabile coniugata a priori per la variabile A, in particolare quando le distribuzioni di A e di B sono membri della famiglia esponenziale. Per ulteriori informazioni cfr. l'articolo sulle distribuzioni composte oppure Liu (1994).[2]

L'implementazione di un campionamento di Gibbs collassato[modifica | modifica sorgente]

Collassamento delle distribuzioni di Dirichlet[modifica | modifica sorgente]

Nei modelli bayesiani gerarchici con variabili categoriali, come il modello di allocazione di Dirichlet latente e vari altri modelli usati nell'elaborazione del linguaggio naturale, è pratica comune "fuori collassare" le distribuzioni di Dirichlet che tipicamente sono usate come distribuzione a priori per le variabili categoriali. Il risultato di questo collasso introduce dipendenze tra tutte le variabili categoriali dipendenti su una data distribuzione a priori di Dirichlet, e la distribuzione congiunta di queste variabili dopo il collasso è una distribuzione di Dirichlet multinomiale. La distribuzione condizionale di una data variabile categoriale in questa distribuzione, condizionata sulle altre variabili, assume una forma estremamente semplice che rende il campionamento di Gibbs anche più facile rispetto al caso in cui il collasso non fosse stato fatto. Le regole sono le seguenti:

  1. Fuori collassare una distribuzione a priori di Dirichlet influisce solo sui nodi genitori e figli della distribuzione a priori. Poiché il nodo generitore è spesso una costante, tipicamente è solo dei nodi figli che ci si deve preoccupare.
  2. Fuori collassare una distribuzione a priori di Dirichlet introduce dipendenze tra tutte le variabili categoriali figlie dipendenti su tale distribuzione a priori, ma non introduce nessuna dipendenza ulteriore tra variabili categoriali figlie. (Questo è importante tenerlo a mente, per esempio, ci sono distribuzioni a priori di Dirichlet multiple collegate tra loro da qualche ditribuzione a priori più generale. Ogni distribuzione a priori di Dirichlet può essere collassata indipendentemente ed influisce solo sulle distribuzioni figlie).
  3. Dopo il collasso, la distribuzione condizionata di una variabile figlia dipendente dalle altre assume una forma semplice: la probabilità di vedere un dato valore è proporzionale alla distribuzione iper a priori per tale valore, e il conteggio di tutti gli altri nodi dipendenti assume lo stesso valore. Nodi non dipendenti dal medesimo valore non devono essere conteggiati. Si noti che la stessa regola si applica in altri metodi di inferenza iterativa, come il metodo variazionale di Bayes o il metodo di massimizzazione del valore atteso; tuttavia, se il metodo implica un mantenimento parziale dei conteggi, allora i conteggi parziali per il valore in questione devono essere sommati per tutti gli altri nodi dipendenti. Talvolta questo conteggio parziale sommato è denominato conteggio atteso (expected count) o con un termine similare. Si noti anche che la probabilità è proporzionale al valore risultante; la probabilità effettiva deve essere determinata mediante normalizzazione su tutti i valori possibili che la variabile categoriale può assumere (cioè sommando il risultato calcolato per ogni possibile valore della variabile categoriale e dividendo per questa somma tutti i risultati calcolati).
  4. Se un dato nodo categoriale possiede dei nodi figli dipendenti (ad esempio quando il nodo è una variabile latente in un modello costituito da una miscela di più distribuzioni), il valore calcolato al passo precedente deve essere moltiplicato per le probabilità condizionali reali (non un valore calcolato proporzionale alla probabilità) di tutti i nodi figli dati i loro nodi genitori. Cfr. l'articolo sulla distribuzione di Dirichlet multinomiale per una discussione dettagliata.
  5. Nel caso in cui l'appartenenza ad un prefissato gruppo dei nodi dipendenti da una distribuzione a priori di Dirichlet possa cambiare dinamicamente in funzione di un'altra variabile (ed esempio una variabile categoriale indicizzata mediante un'altra variabile categoriale latente, come in un modello argomentale (topic model) dove dall'esame di documenti si deve risalire all'argomento in essi trattato), gli stessi conteggi attesi vengono calcolati ancora, ma in maniera accurata in modo da includere l'insieme corretto di variabili. Cfr. l'articolo sulla distribuzione multinomiale di Dirichlet per un'ulteriore discussione, inclusa anche una sul modello argomentale.
Collasso di altre distribuzioni a priori coniugate[modifica | modifica sorgente]

In generale, una distribuzione a priori coniugata può essere fuori collassata, sempre che i nodi figli abbiano distribuzioni coniugate ad essa. I dettagli matematici sono discussi nell'articolo sulle distribuzioni composte. Se è presente un unico nodo figlio, il risultato spesso presume una distribuzione nota. Per esempio, fuori collassando una varianza distribuita come una Gamma Inversa da una rete con un singolo nodo figlio distribuito in forma gaussiana essa darà origine ad una distribuzione di tipo t-Student. (Per quello che importa, fuori collassando entrambe la media e la varianza di un singolo nodo figlio con distribuzione gaussiana darà ancora una distribuzione di tipo t-Student, sempre che le distribuzioni delle prime due siano coniugate, ad esempio gaussiana quella della media, gamma-inversa quella della varianza.)

Se ci sono nodi figli multipli, allora diverranno tutti dipendenti, come nel caso della distribuzione categoriale di Dirichlet. La distribuzione congiunta risultante ha una forma chiusa che assomiglia in qualche modo alla distribuzione composta, nonostante includa il prodotto di un numero di fattori, uno per ciascuno dei nodi figli, in essa.

Inoltre, e maggiormente importante, la distribuzione condizionata risultante dei uno dei nodi figlio dati gli altri (e dati anche i nodi genitore dei nodi collassati, ma non dati i nodi figli dei nodi figli) avrà la stessa densità come la distribuzione predittiva a posteriori di tutti i nodi figli rimanenti. Inoltre, la distribuzione predittiva a posteriori ha la medesima densità della distribuzione composta di base del nodo singolo, nonostante con parametri differenti. La formula generale è data nell'articolo sulle distribuzioni composte.

Per esempio, data una rete di Bayes costituita da un insieme di nodi identici ed indipendenti aventi distribuzione gaussiana con media e varianza aventi distribuzioni a priori coniugate, la distribuzione condizionata di un nodo dati gli altri, dopo aver composto fuori sia la media che la varianza, sarà una distribuzione t-Student. Similmente, il risultato di comporre fuori la distribuzione a priori di tipo gamma di un numero di nodi di tipo poissoniano causa la distribuzione condizionata su un nodo dati gli altri assumere una distribuzione binomiale negativa.

In questi casi dove la composizione produce una distribuzione ben nota, spesso esistono delle procedure efficienti di campionamento, e il loro impiego sarà (ma non necessariamente) più efficiente rispetto alla procedura di collasso, e piuttosto di campionare separatamente le distribuzioni dei nodi figlio e quella a priori. Tuttavia, nel caso in cui la distribuzione composta non sia ben nota, può non essere semplice campionare da essa, poiché in generale non apparterrà alla famiglia delle distribuzioni esponenziali e tipicamente non sarà una funzione logaritmicamente concava (il che renderebbe facile il campionamento mediante il metodo di campionamento con rigetto adattivo (adaptive rejection sampling) e garantirebbe l'esistenza di una forma chiusa).

Nel caso in cui i nodi figli dei nodi collassati abbiano a loro volta dei nodi figli, la distribuzione condizionata di uno di questi nodi figli dati tutti gli altri nodi nel grafo dovrà tenere da conto della distribuzione di questi nodi figli di secondo livello. In particolare, la distribuzione condizionata risultante sarà proporzionale a un prodotto della distribuzione composta come definito sopra, e delle distribuzioni condizionate di tutti i nodi figlio dati i loro nodi genitore (ma non dati i propri nodi figlio). Questo segue dal fatto che la piena distribuzione condizionata è proporzionale alla distribuzione congiunta. Se i nodi figlio dei nodi collassati hanno una distribuzione continua, il campionamento è fattibile, indipendentemente dal fatto che i figli di questi nodi figlio possiedano una distribuzione continua o discreta. In effetti il principio qui implicato è descritto nei dettagli nell'articolo sulla distribuzione multinominale di Dirichlet.

Campionatore di Gibbs con sovrarilassamento ordinato[modifica | modifica sorgente]

  • Un campionatore di Gibbs con sovrarilassamento ordinato (ordered overrelaxation) campiona ad ogni passo un prefissato numero dispari di valori candidati x_j^{(i)} e li ordina, assieme con il singolo valore x_j^{(i-1)} secondo qualche ordinamento ben definito. Se x_j^{(i-1)} è l'elemento s-esimo più piccolo nella lista ordinata, allora x_j^{(i)} è selezionato come l'elemento s-esimo più grande nella lista ordinata. Per maggiori informazioni cfr. Neal (1995).[3]

Altre estensioni[modifica | modifica sorgente]

Il campionamento di Gibbs si può estendere in vari modi. Per esempio, nel caso di variabili casuali le cui distribuzioni condizionate non sono facili da campionare, il campionamento può essere eseguito con il campionamento a strati (slices sampling) oppure l'algoritmo di Metropolis-Hastings. È possibile anche incorporare variabili che non sono variabili casuali ma il cui valore è deterministicamente calcolato da altre variabili. Modelli lineari generalizzati, ad esempio la regressione logistica (ossia modelli di entropia massima), possono essere in tal modo incorporati. (Il software BUGS, per esempio, permette questo tipo di mescolamento dei modelli).

Modalità di insuccesso[modifica | modifica sorgente]

Esistono due modalità in cui il campionamento di Gibbs può avere insuccesso. La prima si verifica quando ci sono isole di stati di elevata probabilità, con nessuno percorso di collegamento tra di loro. Per esempio, consideriamo una distribuzione di probabilità sopra vettori a 2 bit, dove i vettori (0,0) e (1,1) ognuno ha probabilità 1/2, ma gli altri due vettori (0,1) e (1,0) hanno probabilità nulla. Il campionamento di Gibbs rimarrà intrappolato in uno dei due vettotri di probabilità maggiore, non raggiungendo mai l'altro. Più in generale, nel caso di qualsiasi distribuzione definità sopra vettori multidimensionali a valori reali, se due elementi sono perfettamente correlati (o perfettamente anticorrelati), quei due elementi risulteranno indistinguibili per il campionatore di Gibbs.

Un secondo problema può verificarsi quando tutti gli stati hanno probabilità non nulla ed una loro singola piccola "isola" di stati ad elevata probabilità di verificarsi. Per esempio, una distribuzione di probabilità definita sopra un vettore da 100 bit, dove il vettore nullo, cioè quello con tutte le componenti nulle, ha una probabilità pari a 1/2 mentre tutti gli altri vettori sono equiprobabili, ognuno con una probabilità \frac{1}{2(2^{100}-1)}. Se volessimo stimare la probabilità del vettore nullo, sarebbe sufficiente eseguire 100 o 1000 campionamenti dalla distribuzione vera. Otterremmo verosimilmente una risposta vicina a 1/2. Ma probabilmente sarebbero richiesti 2^{100} campioni da un campionamento di Gibbs per ottenere il medesimo risultato. Nessun calcolatore sarebbe in grado di farlo in tempi ragionevoli.

Questo problema si manifesta indipendentemente dal tempo di aggiustamento iniziale dell'algoritmo. Questo accade perché nella distribuzione reale, il vettore nullo occorre per metà del tempo, e queste occorrenze sono casualmente mescolate con le occorrenze dei vettori non nulli. Anche un campione minimo avrà al suo interno il vettore nullo assieme a vettori non nulli. Ma l'algoritmo di campionamento di Gibbs alterna per lunghi periodi la generazione del solo vettore nullo (circa 2^{99} campioni in serie) con la generazione di soli vettori non nulli (in una serie di circa 2^{99} campioni). Perciò la convergenza alla distribuzione vera è estremamente lenta, richiedendo molti di più di 2^{99} passi che non sono dal punto di vista computazionale fattibili in un periodo ragionevole di tempo. La lenta convergenza può essere qui vista come una conseguenza della "maledizione della dimensione", termine coniato da Richard E. Bellman per esprimere un'idea del fenomeno che si manifesta come una progressiva dispersione dei dati all'aumentare della dimensione dello spazio (avente centinaia o migliaia di dimensioni) in cui vengono trattati.

Si noti che un problema come questo può essere risolto campionando in blocco l'intero vettore di 100-bit in una volta. (Questo presume che il vettore da 100-bit sia parte di un insieme più ampio di variabili. Se questo vettore è stato la sola cosa campionata, allora di fatto il campionamento in blocco non equivale al campionamento di Gibbs, il quale per ipotesi sarebbe difficile.)

Software[modifica | modifica sorgente]

Il software OpenBUGS (Bayesian inference Using Gibbs Sampling) fa un'analisi bayesiana dei modelli statistici complessi usando la catena di Markov Monte Carlo.

JAGS (Just another Gibbs sampler) è un programma sotto licenza GPL per l'analisi di modelli gerarchici bayesiani utilizzante la catena di Markov Monte Carlo.

Note[modifica | modifica sorgente]

  1. ^ S. Geman e D. Geman, Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 6, nº 6, 1984, pp. 721–741, DOI:10.1109/TPAMI.1984.4767596.
  2. ^ Jun S. Liu, The Collapsed Gibbs Sampler in Bayesian Computations with Applications to a Gene Regulation Problem in Journal of the American Statistical Association, vol. 89, nº 427, settembre 1994, pp. 958–966, DOI:10.2307/2290921, JSTOR 2290921.
  3. ^ Template:Cite techreport

Riferimenti bibliografici[modifica | modifica sorgente]

Collegamenti esterni[modifica | modifica sorgente]