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