titolo

mercoledì, 19 dicembre 2018

Come costruire un motore Monte Carlo

di mariano tomatis
pubblicato il 9 ottobre 2008

Nel suo Giocati dal caso Nassim Nicholas Taleb racconta di essersi divertito un mondo a risolvere i più astrusi problemi matematici affrontandoli costruendo motori di Montecarlo. Si tratta di un "bricolage" virtuale di gran fascino, di fronte al quale molti matematici storcono il naso (perché i risultati numerici non hanno l'eleganza e la precisione di quelli ottenuti con opportune equazioni), ma possono procurare grande piacere in chi ne mette a punto ogni dettaglio, a ricreare modelli di grande complessità.

Domenica scorsa, all'ingresso di un cinema multisala, ho regalato alla fidanzata un piccolo portachiavi della serie "Pucca"; inserendo 2 euro nel distributore, una pallina saltava fuori rivelando il suo contenuto: uno dei 12 personaggi di una serie dedicata ai segni zodiacali.

Mi sono chiesto quanto si sarebbe dovuto spendere mediamente per ottenere una serie completa di 12 portachiavi, uno per segno: nella migliore delle ipotesi sarebbero bastati 24 euro... ma si sarebbe trattato della situazione ideale, quella in cui non si fosse mai presentato alcun doppione.

Qual è il numero medio di monete che avrei dovuto inserire per ottenere un numero sufficiente di portachiavi a formare una serie completa? (Collect all 12!, recitava un'invitante scritta in alto a destra. Anticipo che la spesa media si aggirerà sui 70 euro).

Ecco come costruire un motore Montecarlo (sviluppato nell'ambiente statistico R) per calcolarlo.

Assegno ad ognuno dei 12 segni zodiacali un numero da 1 a 12. Una serie completa è, quindi, una sequenza di dodici numeri da 1 a 12.

serie <- c(1:12)

Non so con quante serie venga riempito il distributore; supponiamo che vengano introdotte 50 serie. Lo "riempio" virtualmente anch'io nello stesso modo. La funzione rep qui sotto ripete per 50 volte la serie:

distributore <- rep(serie,50)

Diamo ora una bella mescolata al distributore, perché anche nel distributore reale i vari portachiavi si mescolano inesorabilmente:

distributore <- sample(distributore, 50*12, replace = FALSE)

Ho preparato anche 12 scatole, all'inizio tutte vuote, all'interno delle quali collocherò i portachiavi via via acquistati. Ognuna conterrà i portachiavi del segno corrispondente: la prima il segno 1, la seconda il 2, ecc. Solo quando tutte conterranno almeno un portachiavi potrò dirmi soddisfatto e interrompere gli acquisti.

Ognuna delle 12 scatole virtuali conterrà il valore 0 se non ho ancora acquistato alcun segno a lei corrispondente, oppure un numero pari a quanti segni a lei corrispondenti ho acquistato. Se, quindi, la scatoletta 2 conterrà il numero 4 significherà che finora ho acquistato ben quattro portachiavi col segno del toro (il segno 2). Se invece la scatoletta 3 conterrà il numero 0 significherà che non ho ancora trovato il portachiavi dei gemelli.

Inizio mettendo uno zero dentro ognuna delle dodici scatole:

scatola <- rep(0,12)

Posso cominciare ad acquistare portachiavi: devo tenere il conto degli acquisti, e devo andare avanti finché ci sono delle scatole con lo 0 all'interno. In linguaggio logico, proseguo mentre (while) il conteggio delle scatole contenenti lo 0 è positivo - ovvero sum(scatola==0)>0

acquisti<-0
while(sum(scatola==0)>0)
{
   acquisti=acquisti+1
   scatola[distributore[acquisti]]=scatola[distributore[acquisti]]+1
}

Al termine di questo procedimento, in acquisti c'è il numero di portachiavi che ho dovuto acquistare per avere una serie completa.

Il motore è quasi a punto, ma va fatto girare migliaia di volte per produrre un buon risultato. E' ora di innestarlo in una turbina che gli dia almeno diecimila impulsi. Fortunatamente... i soldi spesi per estrarre portachiavi virtuali dal distributore sono anch'essi virtuali!

Predisponiamo un ciclo che lo farà ripetere 10000 volte:

volte <- 10000

Teniamo memoria di quanti portachiavi abbiamo acquistato ad ognuno dei 10000 tentativi. All'inizio fissiamoli tutti a 0, e man mano inseriremo in ognuno il numero di portachiavi acquistati per ottenere una serie completa:

tentativo <- rep(0,volte)

Immergiamo il dispositivo nel ciclo:

for(i in c(1:volte))

Al termine di ogni tentativo assegno al tentativo il numero di acquisti effettuati:

tentativo[i]<-acquisti

E' ora di metterlo in moto: in circa 7 secondi ha completato il suo lavoro. Quanto devo spendere in media per completare una serie?

Chiedo ad R quale sia la media di acquisti nei diecimila tentativi:

mean(tentativo)

Ottengo il numero 36.3213

Poiché ogni portachiavi costa 2 euro, la spesa complessiva è data da:

mean(tentativo)*2

ed è di EUR 72,64 circa.

Chiedendo ad R quale sia stato il tentativo in cui ho speso di meno:

min(tentativo)

scopro che in almeno un'occasione sono stati sufficienti 13 acquisti per completare una serie (per una spesa di 26 euro).

Nell'occasione peggiore, che chiedo così:

max(tentativo)

scopro di aver acquistato 159 portachiavi prima di avere una serie completa, spendendo EUR 318!

Posso chiedere ad R la deviazione standard del numero di acquisti:

sqrt(var(tentativo))

ottenendo circa 12.94

Quest'ultimo valore mi consente di concludere che nel 95% delle occasioni spenderò meno di EUR 123,41.

(mean(tentativo)+1.96*sqrt(var(tentativo)))*2

La distribuzione ottenuta (sull'asse delle ascisse, il numero di acquisti necessari per completare la serie) è la seguente:


hist(tentativo,breaks=100,xlim=c(0,160))

Il motore Montecarlo appena creato può essere scaricato cliccando qui.

Articolo visualizzato 3318 volte

Questo sito non rappresenta una testata giornalistica in quanto non viene aggiornato con cadenza periodica
né è da considerarsi un mezzo di informazione o un prodotto editoriale ai sensi della legge n.62/2001.
Praestigiator è curato da Mariano Tomatis