• View
• Edit
• History
• Print

AltriGeneratoriRigetto

< Integrali con il metodo del rigetto. | Indice | Generatore di numeri casuali gaussiani col metodo di Box e Müller? >

Generatori di numeri casuali: metodo del rigetto

Quando non è possibile calcolare analiticamente l'integrale della distribuzione continua secondo la quale si vogliono estrerre i numeri, oppure non si può invertire analiticamente la relazione che risulta dall'equazione integrale (1) del metodo di trasformazione si può ricorrere al metodo del rigetto.

Questo metodo ricalca il calcolo dell'integrale discusso nella pagina precedente, ed è anaolgo a quello che abbiamo impiegato per ricavare la stima di {$\pi$}. Supponiamo di conoscere la forma analitica della distribuzione {$p(x)$} e di aver determinato sia gli estremi {$(a,b)$} dell'intervallo di {$x$} in cui {$p(x)$} è diversa da zero, sia il massimo {$c$} della funzione {$p(x)$} (per definizione di probabilità quest'ultimo valore è senz'altro minore di 1). Si può allora estrarre due numeri con la distribuzione uniforme standard (compresi tra 0 ed 1), siano {$r_x, r_y$}, e procedere come segue:

• con {$r_x$} si genera {$x=a+(b-a)r_x$}, distribuito uniformemente nell'intervallo {$(a,b)$};
• con {$r_y$} si genera {$y=cr_y$}, distribuito uniformemente nell'intervallo {$(0,c)$};
• il numero {$x$} verrà accettato come estrazione dalla distribuzione {$P(x)$} se {$y \le p(x)$}

Fig. 1 Il metodo del sorteggio con rigetto equivale a valutare l'integrale, ossia l'area {$A$} sotto la curva, con lo stesso metodo che abbiamo visto in precedenza. La probabilità di estrarre un numero compreso tra {$x$} e {$x+dx$} risulterà in questo modo pari al rapporto tra l'area {$dA$} (in blu in figura), sottesa dalla curva {$p(x)$} nell'intervallo {$dx$}, e l'area {$A$} totale. Questo criterio soddisfa la definizione di probabilità per una distribuzione continua, {$dP=p(x)dx$}

Il difetto di questo metodo diventa evidente se il dominio della distribuzione non è limitato, oppure se in un'ampio intervallo il valore della funzione {$p(x)$} è molto minore del minimo, in parole povere se l'area entro il rettangolo rosso, ma sopra la curva blu in Figura 1 è molto maggiore di quella sotto la curva blu. In questo caso il metodo diventa poco efficace perchè si accettano pochi dei valori generati e se ne rigettano molti.

In matlab la generazione avviene con:

r=rand(N,2); % genera N coppie di numeri casuali

x=a+b*r(:,1); % numero casuale uniforme

y=c*r(:,2); % numero casuale uniforme

x=x(find(y<=pdf(x)); % pdf.m contiene la funzione che genera la distribuzione

< Integrali con il metodo del rigetto. | Indice | Generatore di numeri casuali gaussiani col metodo di Box e Müller? >