Supponiamo di voler estrarre dati distribuiti secondo una certa distribuzione con densità di probabilità $f(x)$. In linea di massima, ad esempio col computer, possiamo estrarre campioni dalla distribuzione uniforme; dovremo poi manualmente trasformarli nella distribuzione che ci serve. Un metodo per fare ciò è la cosiddetta estrazione per scarto (rejection sampling in inglese).
L’idea è la seguente. Supponiamo di avere una distribuzione $g(x)$ da cui sappiamo estrarre campioni, e che inoltre soddisfi $f(x) \leq M g(x)$ per una qualche costante $M$ su tutto il dominio di definizione della variabile casuale $x$. Allora estraiamo un campione da $g(x)$ e poi lo accettiamo con probabilità $\frac{f(x)}{Mg(x)}$; se il campione viene rifiutato, ne generiamo un altro sempre da $g(x)$ e lo accettiamo con lo stesso criterio e così via. I campioni così generati avranno la distribuzione di $f(x)$. Un grosso vantaggio di questo metodo è che non è necessario conoscere eventuali costanti di normalizzazione che potrebbero essere difficili da calcolare.
Ora dimostriamo che la tecnica funziona. Prima di tutto dimostriamo che la probabilità di accettare il campione è $1/M$. Infatti accettare con probabilità $p$ è equivalente ad estrarre un numero da una variabile uniformemente distribuita $U$ in $(0,1)$ tale che $U \leq p$. Perciò la probabilità che cerchiamo è
\begin{align*}
P\pqty{U \leq \frac{f(X)}{Mg(X)}} &= \int dx \pqty{U \leq \frac{f(X)}{Mg(X)} \bigg\lvert X=x} P(X=x)=\\
&= \int dx \pqty{U \leq \frac{f(X)}{Mg(X)} \bigg\lvert X=x} g(x)=\\
&= \int dx \frac{f(x)}{Mg(x)} g(x)=\\
&= \frac{1}{M} \int dx f(x) = \frac{1}{M} \\
\end{align*}
Nella prima riga abbiamo usato la legge della probabilità assoluta, poi nella seconda riga il fatto che estriamo $x$ da $g(x)$. Nella terza abbiamo usato il fatto che la probabilità che $U<p$ è proprio $p$, mentre nell’ultima riga l’integrale di $f$ è uno perché $f$ è una densità di probabilità.
Possiamo quindi procedere e dimostrare che i campioni estratti con questo metodo hanno la giusta distribuzione di probabilità. Ovvero
\begin{align*}
P\pqty{X < x \bigg\lvert U \leq \frac{f(X)}{Mg(X)}} &= \frac{P\pqty{ U \leq \frac{f(X)}{Mg(X)} \cap X < x } }{P(U \leq \frac{f(X)}{Mg(X)} )}=\\
&=M \int_{-\infty}^x P\pqty{ U \leq \frac{f(X)}{Mg(X)} \bigg\lvert X=x’ } P(X=x’) dx’=\\
&=M \int_{-\infty}^x \frac{f(x’)}{Mg(x’)} g(x’) dx’=\int_{-\infty}^x f(x’) dx’ \equiv F(x)\\
\end{align*}
dove $F(x)$ è la funzione cumulativa della distribuzione $f(x)$. Nella prima riga abbiamo usato la legge $P(A | B) = P(A \cap B)/P(B)$, nella seconda riga abbiamo esplicitato il denominatore in base al calcolo effettuato prima e abbiamo usato la legge della probabilità assoluta condizionando su $X=x’$; il fatto che $x’ < x$ si manifesta negli estremi di integrazione. Ciò conclude la dimostrazione.
Consideriamo ad esempio $f(x) = A e^{\beta \cos(x)}$ nell’intervallo $0 \leq x < 2\pi$, dove $A$ è una costante di normalizzazione. Questa distribuzione appare ad esempio nella teoria dei campi su reticolo. Non è chiaro come estrarre $x$ con probabilità $f(x)$. Utilizziamo quindi l’estrazione per scarto. Consideriamo la distribuzione uniforme $g(x) = \frac{1}{2\pi}$ su $0 \leq x < 2\pi$ da cui siamo in grado di estrarre campioni. Dobbiamo quindi trovare una costante $M$ tale che $f(x) \leq M g(x)$. In questo caso è abbastanza chiaro che la scelta migliore di $M$ sia semplicemente $M = \mathrm{max} \frac{f(x)}{g(x)} = 2\pi A e^\beta$. Perciò abbiamo $\frac{f(x)}{M g(x)}=e^{\beta (\cos(x)-1)}$. Quindi l’algoritmo per estrarre $x$ generati da $f(x)$ è il seguente: estriamo $x$ uniformemente in $[0,2\pi]$ e poi accettiamo il campione con probabilità $e^{\beta (\cos(x)-1)}$. Notiamo che questo metodo diventa inefficiente per $\beta$ grande: poiché $\cos(x)-1 \leq 0$ allora la probabilità di accettazione del campione decresce esponenzialmente con $\beta$. Notiamo inoltre che non c’è stato bisogno di calcolare esplicitamente la costante di normalizzazione $A$, cosa che è spesso difficile.
Un altro esempio è il seguente. Supponiamo di voler estrarre un punto uniformemente dalla sfera tridimensionale. La densità di probabilità che ci interessa è quindi $f(\vec x)$ dove $f(\vec x) = C$ costante per $\abs{\vec x}\leq 1$ e zero altrimenti. Non è chiaro come estrarre da questa distribuzione. Consideriamo quindi la distribuzione uniforme nel cubo di lato $2$, ovvero $g(\vec x)$ dove $g(\vec x) = D$ costante per $-1 \leq x,y,z \leq 1$ e zero altrimenti (dove $\vec x = (x,y,z)$). Dobbiamo quindi scegliere $M$ tale che $f(\vec x) \leq M g(\vec x)$. Fuori dal cubo dove entrambe sono nulle questa condizione è banale. Tra il cubo e la sfera $f(\vec x)$ è nulla ma non $g(\vec x)$, e quindi la condizione è di nuovo banale. All’interno della sfera, l’ultima situazione da considerare, possiamo quindi scegliere $M$ come $M = C/D$. Quindi abbiamo $\frac{f(\vec x)}{M g(\vec x)} = 1$ dentro la sfera e nulla altrimenti. Perciò per estrarre un punto uniformemente dalla sfera tridimensionale estraiamo tre numeri $x,y,z$ uniformemente dall’intervallo $[-1,1]$ e accettiamo il campione se $x^2+y^2+z^2 \leq 1$.