Campionamento di Gibbs

Voce principale: Catena di Markov Monte Carlo.

In statistica e in fisica statistica, un campionamento di Gibbs o un campionatore di Gibbs è un algoritmo 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); oppure 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), ed è un'alternativa agli algoritmi deterministici impiegati nell'inferenza statistica come l'algoritmo variazionale di Bayes o il metodo della massima verosimiglianza.

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 dopo la morte di Gibbs.^[1]

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^[1] che la sequenza di campioni costituisce una catena di Markov e che la distribuzione stazionaria di tale catena è di fatto la distribuzione congiunta cercata.

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 è 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.

Supponiamo di voler ottenere ${\displaystyle k}$ valori ${\displaystyle x=\{x_{1},\dots ,x_{p}\}}$ da una distribuzione congiunta ${\displaystyle \left.p(x_{1},\dots ,x_{p})\right.}$; si conosce per ciascuna ${\displaystyle x_{j}}$ componente di ${\displaystyle x}$ la probabilità condizionata ${\displaystyle p(x_{j}|x_{1},\dots ,x_{j-1},x_{j+1},\dots ,x_{p})}$.

Denotiamo l'elemento ${\displaystyle i}$-esimo come ${\displaystyle x^{(i)}=\{x_{1}^{(i)},\dots ,x_{p}^{(i)}\}}$ e procediamo come segue:

1. Cominciamo col fissare arbitrariamente il valore ${\displaystyle x_{j}^{(0)}}$ per ogni variabile ${\displaystyle x_{j}}$.
2. Per ogni campione ${\displaystyle i=\{1\dots k\}}$, campioniamo ogni variabile ${\displaystyle x_{j}^{(i)}}$ dalla distribuzione condizionata ${\displaystyle p(x_{j}^{(i)}|x_{1}^{(i)},\dots ,x_{j-1}^{(i)},x_{j+1}^{(i-1)},\dots ,x_{p}^{(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.

Il campione estratto in tal modo approssima la distribuzione congiunta di tutte le variabili. Inoltre, la distribuzione marginale di qualsiasi sottoinsieme di variabili può essere approssimata semplicemente esaminando i valori estratti per quel sottoinsieme di variabili, ignorando le restanti. In aggiunta, il valore atteso di qualsiasi variabile può essere approssimato con una media campionaria.

• 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 sovrarilassamento ordinato; cfr. più sotto.

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

${\displaystyle 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 ${\displaystyle x_{j}}$ e perciò è il medesimo per tutti i valori di ${\displaystyle x_{j}}$. In pratica, per determinare la natura della distribuzione condizionata del fattore ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle 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.

Il campionamento di Gibbs, nella forma dell'algoritmo descritto sopra, definisce una catena di Markov reversibile con la distribuzione invariante desiderata ${\displaystyle p}$. Questo può essere provato nel modo seguente: definiamo ${\displaystyle x\sim _{j}y}$ se ${\displaystyle \left.x_{i}=y_{i}\right.}$ per tutti gli ${\displaystyle i\neq j}$ e sia ${\displaystyle \left.p_{xy}\right.}$ la probabilità di un salto da ${\displaystyle x}$ a ${\displaystyle y}$. Allora, le probabilità di transizione sono

${\displaystyle p_{xy}={\begin{cases}{\frac {1}{d}}{\frac {p(y)}{\sum _{z\sim _{j}x}p(z)}}&x\sim _{j}y\\0&{\text{altrimenti}}\end{cases}}}$

In tal modo

${\displaystyle p(x)p_{xy}={\frac {1}{d}}{\frac {p(x)p(y)}{\sum _{z\sim _{j}x}p(z)}}={\frac {1}{d}}{\frac {p(y)p(x)}{\sum _{z\sim _{j}y}p(z)}}=p(y)p_{yx}}$

in quanto ${\displaystyle 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 ${\displaystyle \left.g\right.}$.

In pratica, il suffisso ${\displaystyle \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).

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.

• Un campionamaneto di Gibbs a blocchi raggruppa assieme due o più variabili e campiona la loro distribuzione congiunta condizionata su tutte le variabili restanti, piuttosto che campionare da ogni singola variabile individualmente. Per esempio, in un modello di Makov nascosto, un campionamento di Gibbs a blocchi può campionare da tutte le variabili latenti costruendo la catena di Markov in un solo colpo mediante l'algoritmo forward-backward.

• 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]

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 distribuzione 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.

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.

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

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).

Può succedere che due o più variabili xi da cui campionare siano fortemente correlate (positivamente o negativamente), se questo succede, al momento di campionare da ciascuna di queste variabili, il condizionamento alle altre con cui è correlata farà sì che il valore campionato tenda a non spostarsi mai dalla regione in cui si trova. Nel caso limite di due variabili perfettamente correlate, al momento di campionarne una, dato che la sua distribuzione si considera condizionata all'altra variabile, sarà estratto un valore perfettamente uguale a quello dell'altra variabile, e viceversa, cosicché qualunque sia il valore iniziale della seconda di queste variabili, tale valore non cambierà mai per tutta la catena, impedendone la convergenza. Si noti che se si conosce (o si può simulare) la distribuzione condizionata di tutto il blocco delle variabili tra loro correlate, il problema si può facilmente risolvere impostando un Gibbs sampler a blocchi.

Un secondo problema è conosciuto come "maledizione della dimensionalità" (curse of dimensionality, termine coniato da Richard Bellman e usato in una varietà di contesti diversi ma simili) e si verifica per problemi in cui la dimensionalità delle x è molto alta. In questi casi la regione a maggior densità tende a perdersi all'interno dello spazio delle x e il Gibbs sampler risulta quindi fortemente autocorrelato. Per fare un esempio semplice, una distribuzione di probabilità è definita sopra un vettore da 100 bit (100 variabili dicotomiche), 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à ${\displaystyle (2(2^{100}-1))^{-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 ${\displaystyle 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. 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 ${\displaystyle 2^{99}}$ campioni in serie) con la generazione di soli vettori non nulli (in una serie di circa ${\displaystyle 2^{99}}$ campioni). Perciò la convergenza alla distribuzione vera è estremamente lenta, richiedendo molto più di ${\displaystyle 2^{99}}$ passi, che non sono dal punto di vista computazionale fattibili in un periodo ragionevole di tempo.

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.

• George Casella e Edward I. George, Explaining the Gibbs sampler, in The American Statistician, vol. 46, n. 3, 1992, pp. 167–174, DOI:10.2307/2685208, JSTOR 2685208. (Contiene un sommario riassuntivo e molti riferimenti.)
• Alan E. Gelfand e Adrian F. M. Smith, Sampling-Based Approaches to Calculating Marginal Densities, in Journal of the American Statistical Association, vol. 85, n. 410, 1990, pp. 398–409, DOI:10.2307/2289776, JSTOR 2289776, MR 1141740.
• Andrew Gelman, John B. Carlin, Hal S. Stern, and Donald B. Rubin. Bayesian Data Analysis. London: Chapman and Hall. First edition, 1995. (See Chapter 11.)
• C.P. Robert and G. Casella. "Monte Carlo Statistical Methods" (second edition). New York: Springer-Verlag, 2004.
• David A. Levin, Yuval Peres, and Elizabeth L. Wilmer. "Markov Chains and Mixing Times". https://web.archive.org/web/20100516063638/http://www.uoregon.edu/~dlevin/MARKOV/
• Bolstad, William M. (2010) Understanding Computational Bayesian Statistics, John Wiley ISBN 978-0-470-04609-8
• Christopher M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006, ISBN 0-387-31073-8.

• The OpenBUGS Project - Bayesian inference Using Gibbs Sampling, su openbugs.info.
• Una semplice spiegazione del campionamento di Gibbs può essere trovata nelle pp. 370–371 del volume del professor MacKay's Information Theory, Inference, and Learning Algorithms, disponibile per libera consultazione al sito inference o al sito dell'università di Toronto.
• Un'applicazione pratica del campionamento di Gibbs in genomica., su bayesweb.wadsworth.org. URL consultato il 27 ottobre 2011 (archiviato dall'url originale il 4 novembre 2011).

Portale Statistica: accedi alle voci di Wikipedia che trattano di statistica