27 febbraio 2009

Numeri pseudo-casuali e software statistici

La generazione di numeri pseudo-casuali è un argomento di fondamentale importanza per le applicazioni in statistica. La "bontà" dell'algoritmo utilizzato in un software è quindi decisivo per la validità delle simulazioni eseguite mediante esso. In questo breve post intendo soffermermi sugli algoritmi impementati in R, SAS ed SPSS.
L'argomento è stato trattato ampiamente dagli studiosi e continua ad attirare l'attenzione di statistici ed informatici (...ammesso che in tale ambito sia possibile fare questa distinzione tra le due discipline...), ed ovviamente gli sviluppatori dei software sopra citati non possono che seguire gli sviluppi della scienza.
I numeri peseudo-casuali sono chiamati in tale modo appunto perché sono "deterministici" ma "sembrano" casuali. Sono generati da un algoritmo implementato in un calcolatore (quindi deterministisco) ma soddisfano una serie di proprietà che li fanno "sembrare" casuali per gli scopi prefissati, quindi per la simulazione di un qualche fenomeno. L'algoritmo più famoso è il Generatore Lineare Congruente (Linear Congruential Generators) riassunto dalla seguente funzione di trasferimento:
f(x) = (ax + c) mod m .
mod è la funzione che fornisce il resto della divisione intera, ossia:
x mod y = y ( x/ y - [x/y] )
dove [h] è la parte intera inferiore di h. Quindi ad esempio 11 mod 3 = 3 * (11/3 - [11/3]) = 3*(3.6667 - 3) = 2.
Qundi si fissano a, c ed m, si parte da un valore inziale X_0 e si ottiene X_1 = f(X_0) . A questo punto si genererà il numero pseudo-casuale (in [0, 1[ ) mediante una funzione output (la più semplice è X_j /m). La funzione di trasferimento più semplice e più conosciuta è invece quella di Park-Miller, in cui si impone semplicemente c=0.
Per una trattazione teorica molto chiara ed esaustiva si legga qui.
La tipologia di algoritmo appena citato è stato nel tempo criticato e praticamente sostituito dal generatore di Mersenne-Twister (M-T).
Ovviamente, poiché R è sempre "avanti" :-) , quest'ultimo è l'algoritmo di default già da tempo, mentre gli altri software sono stati aggiornati da un po'.
In SAS infatti, fino alla versione 9.0, l'unico algoritmo disponibile era il Generatore Lineare Congruente nella versione di Park-Miller con i seguenti parametri:
c=0
a=397204094
m=2^31-1 .
Solo dalla 9.1 è finalmente disponibile anche la funzione RAND che appunto si basa sull'algoritmo di M-T!!!! (chi vuole approfondire può leggere qui).
Mi dilungo sul SAS perché trovo molto macchinoso il procedimento di gestione delle sequenze dei numeri casuali. Utilizzando la funzione RANUNI(seed) in un datastep, cambiare il seed (ossia il valore di x_0 nel generatore) non avrà influenza sulla sequenza di numeri casuali poiché l'estrazione avverrà sempre dallo stesso stream (ossia flusso, ma in inglese fa più figo :-) ). Per avere un controllo della sequenza sarà necessario ricorrere a CALL RUNUNI(seed) dove il seed deve essere una variabile e non una costante.
Si provi ad eseguire il seguente codice per capire la distinzione:

data test;
do i=1 to 10;
x = RANUNI(1);
output;
end;

run;

data test_i;
do i=1 to 10;
x = RANUNI(i);
output;
end;
run;


data test_call;
do i=1 to 1;
call RANUNI(i,x); /*con ranuni(1,x) va in errore*/
output;
end;
run;


Per quanto riguarda infine l'SPSS, si faccia riferimento al menu " Trasforma / Generatori numeri casuali " per scoprire il tipo impostato. Consiglio ovviamente di utilizzare il generatore di Mersenne-Twister con inizializzazione casuale (ossia derivata dall'orologio di sistema). Più semplice forse farlo direttamente aprendo un file di sintassi .sps e settando mediante codice:
SET RNG=MT MTINDEX=RANDOM.

Nessun commento: