Simulazione del paradosso di S. Pietroburgo

[albert]

Il grafico riportato nell’articolo http://www.riflessioni.it/scienze/pa...ietroburgo.htm è stato ottenuto con una semplice simulazione al computer. Si è utilizzata la funzione pseudocasuale Unix “rand”, l’uscita di un numero dispari corrispondeva a testa.

main(){

unsigned int giocate, num;
int testa = 1;
float media,total=0;

srand(31415);

for(giocate=1;giocate<20000; ){
num=rand();
if(num%2) testa <<= 1;
else {
total += testa;
media = total/giocate;
printf("%d; %5.3f\n", testa, media);
testa=1;
giocate+=1;
}
}
}

L’output è stato ridiretto in un file poi aperto con excel, con cui si è creato il grafico. Il valore finale è di 7.4, in accordo con il valore atteso di circa 7.5, secondo la formula:

E=½ + ½ log2 n

dove n è 20000.

Estendendo la simulazione ad un numero molto più alto di giocate (100 milioni) ci si aspetterebbe un valore di circa 14. Ho però notato che, con questa funzione pseudocasuale, il valore non saliva oltre 10.3. Credo che il problema sia la insufficiente bontà della funzione, che non permette più di 16 occorrenze di numeri pari (o dispari). Esaminando infatti le vincite di 128 (7 occorrenze) o più nel file delle 100 milioni di giocate, si osserva il comportamento riportato nella tavola (seconda colonna numero effettivo di presenze, terza colonna numero previsto):

128 382737 390625
256 190643 195312
512 104168 97656
1024 43297 48828
2048 23456 24414
4096 7803 12207
8192 8373 6103
16384 2518 3051
32768 7916 1525
65536 536 762

A partire da 4096 (12 occorrenze) si cominciano a notare degli scostamenti abbastanza rilevanti. Non si osservano valori pari a 131072 (17 occorrenze) o superiori.

Qualcuno ha qualche idea su come realizzare una funzione pseudocasuale migliore??

[epicurus]

Originalmente inviato da albert

Si è utilizzata la funzione pseudocasuale Unix “rand”, l’uscita di un numero dispari corrispondeva a testa.

main(){

unsigned int giocate, num;
int testa = 1;
float media,total=0;

srand(31415);

for(giocate=1;giocate<20000; ){
num=rand();
if(num%2) testa <<= 1;
else {
total += testa;
media = total/giocate;
printf("%d; %5.3f\n", testa, media);
testa=1;
giocate+=1;
} } }

[...]

Qualcuno ha qualche idea su come realizzare una funzione pseudocasuale migliore??

Ciao Alberto. Innanzitutto dovresti sostituire srand(31415) con srand(time(NULL)), altrimenti ogni volta avrai lo stesso seme, quindi la stessa sequenza pseudocasuale.

Per quanto riguarda a che funzione pseudocasuale migliore consigliarti, bisognerebbe prima sapere che funzione sia effettivamente la tua rand, cioè che algoritmo implementa. Comunque prova a dare un occhio all'algoritmo Mersenne twister, o a Blum Blum Shub.

[albert]

Originalmente inviato da epicurus

Ciao Alberto. Innanzitutto dovresti sostituire srand(31415) con srand(time(NULL)), altrimenti ogni volta avrai lo stesso seme, quindi la stessa sequenza pseudocasuale.

Sì, hai ragione, ma avevo un motivo ... in questo modo ho potuto verificare se librerie differenti davano lo stesso risultato, e quindi usavano lo stesso algoritmo. Su HP-UX e sulle librerie "C" della Zortech che uso sul mio PC da' esattamente la stessa sequenza. Non è però documentato l'algoritmo. Su Linux RedHat (ho provato adesso) dà una incredibile vincita di 524288 dopo una decina di giocate, per cui i valori dei run "corti" sono molto più alti del previsto. A 100 milioni dà 13.6, abbastanza concorde con la teoria. Ma visto la clamorosa vincita iniziale (che sbancherebbe subito il casinò) non saprei dire se sia un "buon algoritmo"

Originalmente inviato da epicurus

Per quanto riguarda a che funzione pseudocasuale migliore consigliarti, bisognerebbe prima sapere che funzione sia effettivamente la tua rand, cioè che algoritmo implementa. Comunque prova a dare un occhio all'algoritmo Mersenne twister, o a Blum Blum Shub.

purtroppo nessuna libreria documenta l'algoritmo, sigh