Probabilit a e computer
Modello probabilistico e spesso un compromesso tra
esigenza di aderenza a fenomeno reale e trattabilita
matematica ⇒ modello troppo semplificato ma che puo
essere analizzato matematicamente
Evoluzione dei mezzi di calcolo rende possibile modello
piu accurato, da analizzare attraverso uno studio di
simulazione
Simulazione rende anche possibile studiare fenomeni
virtuali o organizzare esperimenti virtuali
Esempio: Coda ad uno sportello
Si prevede di aprire uno sportello che dia informazioni tra
le 9 e le 19 ogni giorno.
Ci si aspetta che ogni giorno si presentino allo sportello
circa 50 persone (clienti).
Si prevede che rispondere ad un cliente richieda un
tempo casuale con media 10 e deviazione standard 2.
Si pensa di non accettare piu clienti dopo le 19 ma di
servire comunque quelli in attesa.
Tipiche domande alle quali si vuole ottenere una risposta:
1
1. tempo medio di apertura effettiva dello sportello?
2. frazione di giorni in cui lo sportello ancora aperto alle
19,30?
3. tempo medio di attesa per un cliente?
4. quanti clienti in media serviti ogni 30 minuti?
5. effetto di respingere i clienti quando la coda e piu
lunga di 5 clienti?
Costruzione di un modello probabilistico “utile” richiede
ipotesi ragionevolmente accurate.
• modello del processo degli arrivi: tasso di arrivo
costante? Variabile con le ore del giorno (piu
realistico)
• distribuzione del tempo di servizio di un cliente?
Diverse distribuzioni in dipendenza del numero di
clienti in attesa o delle ore del giorno?
• dipendenza dai giorni della settimana?
2
Modello risultante puo essere troppo complicato da
trattare per via analitica.
Simulazione programma comportamento probabilistico
del sistema su un computer utilizzando “numeri casuali”
per realizzare i vari effetti probabilistici (valori delle
variabili aleatorie in gioco) e registrando esito di tali effetti
su aspetti del sistema complessivo di interesse.
Computer permette di replicare “esperimento virtuale”
un gran numero di volte e ottenere stime delle quantita
interessanti tramite legge grandi numeri.
3
Generazione di numeri casuali
Lancio moneta, dado, esito della roulette, estrazione di
carte da un mazzo: sistemi “meccanici” per simulare
realizzazioni di variabili aleatorie discrete
Lancio moneta (equilibrata) puo generare numeri casuali
con arbitrario numero di cifre binarie, approssimazione di
numeri “uniformi” in (0,1) (realizzazioni di v.a. uniforme)
Poco pratico
Con computer numeri “pseudo-casuali”: generati da
algoritmo deterministico, ma che “sembrano” realizzazioni
di variabili aleatorie uniformi e indipendenti
Approccio piu comune: congruenziale moltiplicativo
x0 valore intero iniziale (seme)
xn ≡ axn−1 mod m
a e m interi positivi fissati
xn ∈ {0, 1, . . .m − 1}
xn/m valore approssimato (?) di un numero estratto da
U(0, 1)
4
a e m scelti in modo che:
1. per ogni “ragionevole” scelta di x0,{
x1
m , . . . , xn
m
}
si
comporta come una successione di valori osservati di
n v.a. indipenendenti e U(0, 1)
2. il “periodo” della successione e grande (la succes-
sione prima o poi si ripete, sicuramente dopo m
iterazioni)
3. il calcolo e veloce
Per un computer con parole di 32 bit una buona scelta e
m = 231 − 1 = 2147483647 e a = 75 = 16807
Proprieta statistiche della successione valutate attraverso
una serie di test
231 − 1 numero primo (scelta ... non casuale!)
periodo pieno = m sufficientemente lungo in pratica
Generatore congruenziale misto
xn ≡ (axn−1 + c) mod m
Ogni sistema operativo ha in librerie di sistema procedura
per generare numeri pseudo-casuali
5
Calcolo di un integrale definito tramite numeri pseudo-
casuali
Una delle prime applicazioni dei numeri pseudo-casuali:
metodo Montecarlo (all’inglese Monte Carlo)
I =
∫ 1
0
g(x)dx
Se U ∼ U(0, 1) ⇒ I = E[g(U)]
Se U1, U2, . . . , Un i.i.d ∼ U(0, 1), da legge forte dei
grandi numeri:
1
k
k∑
i=1
g(Ui) → E[g(U)] = I per k → ∞
Generando un numero elevato di numeri pseudo-
casuali u1, u2, . . . , un, la media aritmetica dei g(ui)
approssima I
6
> # Esempio di calcolo di integrale definito con metodoMontecarlo
> with(stats):
> with(random):
> g:=proc(x) exp(−sin(x)) end proc: #funzione integranda
> int(g(x),x=0..1); #integrazione definita
0
1
e−sin x( )( )
dx
>
> evalf(%); #valuta numericamente
0.6513249125
> media:=proc(g,n)
local i, m:
m:=0:
for i from 1 to n do
m:=m+g(uniform[0,1](1)):
od:
m:=m/n:
end proc:
> media(g,10);
0.6511945812
> media(g,100);
0.6527562463
> media(g,1000);0.6467233088
> media(g,10000);
0.6492964283
> media(g,100000);
0.6511191300
Deviazione standard della media aritmetica ≍ 1√
n, errore
nella stima ≍ 1√
n
7
Se integrale non in [0,1], si puo operare cambio di
variabile
Metodo Montecarlo poco utile per integrali semplici, piu
utile per integrali multipli
I =
∫ 1
0
∫ 1
0
· · ·
∫ 1
0
g(x1, x2, . . . , xn)dx1dx2 · · · dxn
I = E[g(U1, U2, . . . , Un)] ≃k∑
j=1
g(uj1, . . . , u
jn)/k
Ui i.i.d. U(0, 1) e uji numeri pseudo-casuali
Es. Stima Montecarlo di π
(X, Y ) punto casuale uniforme sul quadrato di vertici
(±1,±1)
Probabilita p che (X, Y ) ∈ cerchio iscritto di raggio 1 ?
p = Area del cerchioArea del quadrato
= π4
Se X e Y sono indipendenti e U(−1, 1), f(x, y) =
f(x)f(y) = 12
12 = 1
4 , x, y ∈ (−1, 1) ⇒ (X, Y )
uniforme sul quadrato
Se U ∼ U(0, 1) ⇒ 2U − 1 ∼ U(−1, 1)
8
Da coppie di numeri pseudo-casuali u1, u2 ottengo punti
pseudo-casuali (2u1 − 1, 2u2 − 1) nel quadrato
I(x, y) = 1 se x2 + y2 ≤ 1, 0 altrimenti
⇒ p = E(I) = P (X2 + Y 2 ≤ 1)
⇒ p ≃ frazione di punti pseudo-casuali che cade nel
cerchio
9
> # Esempio di calcolo di pigreco con metodo Montecarlo
> with(stats):
> with(random):
> evalf(Pi); #valuta numericamente3.141592654
> g:=proc(x,y) if (x^2+y^2 <= 1) then 1. else 0 end if: end proc:
> media:=proc(g,n)
local i, m , u :
m:=0:
for i from 1 to n do
u:=uniform[0,1](2):
m:=m+g(2*u[1]−1,2*u[2]−1):
od:
m:=4*m/n:
end proc:
> media(g,10);
2.800000000
> media(g,100);
3.000000000
> media(g,1000);3.168000000
> media(g,10000);
3.151600000
> media(g,100000);
3.142200000
10
Generazione di numeri pseudo-casuali non uniformi
Da numeri pseudo-casuali uniformi simulazione variabili
aleatorie con altre distribuzioni
Simulazione di variabili aleatorie discrete
P (X = xj) = pj , j = 1, 2, . . . ,∑
j
pj = 1
Se U ∼ U(0, 1)
P (X = xj) = P
(
j−1∑
i=1
pi ≤ U <
j∑
i=1
pi
)
= pj
⇒ a u numero uniforme associo xj se
j−1∑
i=1
pi ≤ u <
j∑
i=1
pi
algoritmo:
genera u
se u < p1 stop con x1
se u < p1 + p2 stop con x2...
11
Se xj = j, j = 1, . . . , n e pj = 1n (v.a. uniforme
discreta) l’algoritmo si ferma quando j−1n ≤ u < j
n
restituendo j, cioe [nu] + 1 ([x] parte intera di x) -
calcolo immediato
Simulazione del lancio di una moneta
p probabilita di una “testa” , se u < p testa, altrimenti
croce
Simulando n lanci di una moneta e sommando i
“successi” simulo v.a. binomiale B(n, p)
Non molto efficiente
Oltre a proprieta statistiche, importante efficienza degli
algoritmi
12
algoritmo piu efficiente per binomiale:
P (X = i+1) =n!
(i + 1)!(n − (i + 1))!p(i+1)(1−p)n−(i+1)
=n − i
i + 1
n!
i!(n − i)!
p
1 − pp
i(1−p)n−i =n − i
i + 1
p
1 − pP (X = i)
Applicazione algoritmo generale:
Passo 1: genera u
Passo 2: c = p/(1−p), i = 0, P = (1−p)n, F = P
Passo 3: se u < F stop con i
Passo 4: P = [c(n − i)/(i + 1)]P, F = F + P, i =
i + 1
Passo 5: Vai a passo 3
In media 1+np ricerche; se p > 1/2, conviene simulare
B(n, 1 − p) e sottrarre da n
13
Simulazione della rovina del giocatore
> # Simulazione rovina del giocatore
> with(stats):
> with(random):
> vincita:=proc(p) #1 o −1 con prob p e 1−p if ( uniform[0,1](1) < p) then 1 else −1 end if:end proc:
> gioca:=proc(m,n,p) #calcola se rovina o vittoria e duratagioco local rovina, durata, capitale: durata:=0: capitale:=m: while ( capitale <> 0 and capitale <> m+n ) do capitale:= capitale + vincita(p): durata:= durata + 1: end do: if ( capitale = 0 ) then rovina:=1 else rovina:=0 end if: [rovina,durata]end proc:
> medie:=proc(m,n,p,N) # calcola frequenza rovina e duratamedia # su N ripetizioni
local i, L, prob, esito:
prob:=0: L:=0:
for i from 1 to N do esito:= gioca(m,n,p):
prob:=prob+esito[1]: L:=L+esito[2]:
od:
prob:=prob/N: L:=L/N: evalf([prob,L]):
end proc:
14
> medie(10,10,0.5,10000); # gioco equo con pari capitaleiniziale
0.4984000000, 100.6732000[ ]
> # da confrontare con n/(m+n) = 1/2
# e m*n =10*10 = 100
> medie(1,100,0.4,10000); #ubriaco ad 1 passo dal burrone ea 100 # passi da casa attratto dalburrone
1., 4.875200000[ ]
> subs(m=1,n=100,r=0.6/0.4,( r^m − r^(m+n) )/(1 − r^(m+n) ) );
# probabilitàeffettiva
1.000000000
> subs(m=1,n=100,q=0.6, p=0.4, (m+n)/(p−q) * (1−(q/p)^m) /(1−(q/p)^(m+n)) − m/(p−q) ) # durata mediaeffettiva
5.000000000
> medie(1,100,0.5,10000); # idem ma senza attrazione
0.9905000000, 96.80590000[ ]
> # da confrontare con100/101=0.9901 # e 1*100=100
> medie(1,100,0.6,10000); #idem ma con attrazione verso casa0.6742000000, 159.5712000[ ]
> subs(m=1,n=100,r=0.4/0.6,( r^m − r^(m+n) )/(1 − r^(m+n) ) ); # probabilitàeffettiva
0.6666666667
> subs(m=1,n=100,q=0.4, p=0.6, (m+n)/(p−q) * (1−(q/p)^m) /(1−(q/p)^(m+n)) − m/(p−q) ); #durata mediaeffettiva
163.3333333
15
Simulazione di variabili aleatorie continue
Se U ∼ U(0, 1) e F funzione di ripartizione continua
⇒
X = F−1(U) ∼ F
Algoritmo della trasformazione inversa per simulare una
v.a. con distribuzione F
Genera u
Calcola F−1(u)
esempio
Simulare X ∼ Exp(1)
F (x) = 1 − e−x ⇒ x = F−1(u) = −log(1 − u)
numero pseudo-casuale esponenziale
se U ∼ U(0, 1) ⇒ 1 − U ∼ U(0, 1), quindi basta
porre x = −log(u)
Per simulare X ∼ Exp(λ) x = − 1λ log(u)
Simulazione di un processo di arrivi con intervallo tra un
arrivo e l’altro descritto da v.a. esponenziale
Metodo di accettazione-rifiuto
Se possibile simulare in modo efficiente una v.a. con
16
funzione di densita g(x), possiamo simulare qualsiasi
altra v.a. dotata di densita f(x) se ∃c tale che
f(x)g(x) ≤ c ∀x
Passo 1: Genera Y con densita g
Passo 2: Genera U uniforme
Passo 3: Se U ≤ f(Y )cg(Y ) stop con Y altrimenti vai al
Passo 1
Si dimostra che la variabile prodotta dall’algoritmo ha
densita f e che il numero di iterazioni richieste e una v.a.
geometrica con media c
Esempio
Simulare una v.a. Normale standard
se X ∼ N(0, 1), Z = |X | ha densita f(z) =2
√
2πe−z2/2, 0 < z < ∞
g(z) = e−z densita esponenziale con media 1
f(z)g(z) =
√
2/πez−z2/2
massimo in 1 ⇒ c = f(1)/g(1) =√
2e/π
⇒ f(z)cg(z) = exp
{
− (z−1)2
2
}
17
simulazione di Z :
Passo 1: Genera Y ∼ exp(1)
Passo 2: Genera U ∼ U(0, 1)
Passo 3: Se U ≤ exp{−(Y − 1)2/2} stop con Y
altrimenti vai a Passo 1
X ∼ N(0, 1) puo essere simulata scegliendo Z o −Z
con probabilita 1/2.
18
Sistemi di funzioni iterate
I sistemi di funzioni iterate costituiscono un metodo
per costruire immagini frattali (M. Barnsley, Fractals
Everywhere, 1988).
I frattali sono composti dall’unione di copie di se stessi
ottenute tramite la trasformazione attraverso opportune
funzioni.
L’esempio piu famoso e il triangolo di Sierpinski
• sia dato un triangolo equilatero con base parallela
all’asse orizzontale (prima immagine)
• si riduca ogni lato del triangolo di 1/2, si facciano
due copie e si posizionino i tre triangoli come nella
seconda immagine
• si ripeta il passo precedente per ciascuno dei triangoli
(immagine 3 e seguenti)
Immagine “autosimile”, fatta dall’unione di tre copie
dell’immagine traformate come sopra
Costruzione tramite IFS
19
x0 = 0, y0 = 0
i punti successivi sono ottenuti applicando a caso (con
uguale probabilita) una delle tre seguenti trasformazioni
(affini)
xn+1 = 0.5xn
yn+1 = 0.5yn
xn+1 = 0.5xn + 0.5
yn+1 = 0.5yn + 0.5
xn+1 = 0.5xn + 1
yn+1 = 0.5yn
Le tre trasformazioni producono i punti rispettivamente
evidenziati in giallo, rosso, verde nella figura
20
Catena di Markov (con spazio degli stati continuo)
Le funzioni di trsaformazione non necessariamente
affini, ma devono essere contrazioni, cioe rendere i punti
trasformati piu vicini.
La forma del frattale e di conseguenza fatta di un certo
numero di copie (eventualmente con sovrapposizioni)
ridotte del frattale e cosı ogni copia e fatta di sue copie
ridotte e cosı via.
S = ∪Ni=1fi(S)
I sistemi IFS sono usati per la compressione di immagini.
Descrizione con pochi parametri.
Problema difficile trovare in generale il sistema che
determina una data immagine (brevetti di Barnsley)
IFS per immagine di foglia di felce:
21
La figura e formata dall’unione di 4 copie di se stessa
(compreso lo stelo - verde)
22
Le 4 trasformazioni corrispondenti sono:
xn+1 = 0
yn+1 = 0.16yn stelo - verde
xn+1 = 0.2xn − 0.26yn
yn+1 = 0.23xn + 0.22yn + 1.6 rosso
xn+1 = −0.15xn + 0.28yn
yn+1 = 0.26xn + 0.24yn + 0.44 blu
xn+1 = 0.85xn + 0.04yn
yn+1 = −0.04xn + 0.85yn + 1.6 azzurro
Le probabilita delle trasformazioni sono rispettivamente
0.01, 0.07, 0.07, 0.85.
x0 = y0 = 0
Le probabilita non influenzano la forma, ma solo la
frequenza dei punti nelle varie regioni
23