Numeri casuali -...

Generazione di numeri casuali e pseudo-casuali 1

Numeri casuali

Nei problemi di Simulazione si presenta spesso la necessità di generare determinazioni, ovvero numeri, che riproducano un campione estratto da una certa distribuzione.

A tale scopo risulta quindi necessario disporre di un sistema che generi numeri casuali.

Vedremo come sia sufficiente saper generare successioni di numeri casuali uniformemente distribuiti nell’intervallo (0,1) eforniremo alcuni metodi per poter riprodurre, da questi, i valori di qualsiasi altra v.a. (variabile aleatoria).

Cosa si intende per “successione di numeri casuali” ?Due sono le definizioni disponibili:• “classica”• “moderna”


Una sequenza di numeri casuali è costituita da numeri interi appartenenti all’insieme {1,…,N} e generati tramite ripetizione di un esperimento casuale(prove fra loro indipendenti ed equiprobabili)

q Definizione “classica” ( )Rett 1,u X N→ ∼

Una sequenza di numeri casuali è costuita da numeri reali appartenenti all’intervallo (0,1) indipendenti edequiprobabili

q Definizione “moderna” ( )Unif 0,1u X→ ∼

Numeri casuali

In entrambi i casi, le sequenze di n.c. non devono esserenè prevedibili nè riproducibili


Generatori di numeri casuali

Per generare manualmente successioni di numeri casuali “classici” venivano e vengono usati dispositivi meccanici:

Esempio:

urna contenente palle numerate

Estrazioni a casoripetute

con o senzareinserimento

Successione di numeri casualisecondo la def. classica

Altri esempi: • Dadi “regolari”• Monete “bilanciate”• Mazzi di carte “ben mescolati”

1

2

3

..N


Generazione di numeri casuali

Un altro metodo per generare successioni di numeri casuali consiste nell'uso di dispositivi naturali ovvero nella registrazione di un qualche processo fisico che la teoria postula avvenire in modo casuale.

Ad esempio tramite la codifica del rumore termico di un'apparecchiatura elettronica è possibile costruire sequenze dinumeri casuali di pronto utilizzo o memorizzabili in tavole di numeri casuali per un utilizzo posteriore.

In entrambi i casi il metodo dei dispositivi naturali è ormai completamente abbandonato per la sua complessità d'uso con il calcolatore.Le tavole precompilate, infatti, richiedono molta memoria per il loro immagazzinamento mentre l'interfacciamento diretto dei dispositivi rallenta significativamente la velocità di esecuzione dei programmi.


Numeri pseudocasuali

L'unico modo per eliminare tutte queste difficoltà è quello di far generare i numeri casuali direttamente al calcolatoremediante appositi algoritmi matematici.

Le sequenze di numeri così prodotti vengono denominate pseudocasuali in quanto essendo prodotte da algoritmi deterministici esse, al contrario dei numeri "effettivamente" casuali, non solo sono esattamente prevedibili ma possono anche essere riprodotte in maniera identica quante volte si vuole semplicemente ripetendo il procedimento di calcolo.

Il loro utilizzo viene però pienamente giustificato dal fatto che non sono necessari numeri che soddisfino tutti i criteri formali di casualità ma numeri che abbiano le stesse proprietà di quelli casuali: poco importa come siano stati generati, basta che, se sottoposti a verifica, risultino indistinguibili da quelli veramente casuali in termini di buon adattamento e indipendenza.


Test di casualità

In letteratura esiste una vasta gamma di test di casualità:

Ø “generali”: bontà di adattamento a una f.r. qualunque

Ø test Chi Quadrato

Ø test di Kolmogorov-Smirnov

Ø “specifici”: bontà di adattamento alla f.r. dell’

Ø test di equidistribuzione

Ø “tipici”

Ø Up and Down test (Run test, Wald-Wolfowitz, 1940)

Ø Gap test (per le sequenze di numeri interi)

( )1,0Unif


Generatori di numeri pseudocasuali

Vi sono vari metodi per generare numeri pseudocasuali:

Ø Middle-Square Method (J. von Neumann, 1945/50)

Ø Metodi Congruenziali Lineari (Lehmer, 1951)

Ø Altri metodi (Gentle, 1998)

Le sequenze ottenute con tali metodi hanno:

distribuzione uniforme

In particolare, i metodi congruenziali forniscono sequenze di numeri pseudocasuali con:

distribuzione uniforme "continua" in (0,1)



q Metodo Middle-Square

Ø Si vuole produrre un numero pseudo-casuale di N cifre

Ø Si assegna ad arbitrio un numero X di N cifre (Seme)

Ø Si considerano le N cifre centrali di 2X

Esempio:

61,2 == XN

( ) 7237212 ⇒=X

( ) 18518472 2 ⇒=

( ) 2432418 2 ⇒=

( ) 7657624 2 ⇒=

…



Il metodo Middle-Square non è un buon generatore di numeri pseudocasuali

Ø Sensibile alla scelta del seme

Esempio:

60,2 == XN ( ) 6036002 ⇒=X sequenza costante

33,2 == XN ( ) 0810892 ⇒=X non è di N cifre

Ø Qualora fornisca 0 (00, 000, ….) termina la sequenza



q Metodi congruenziali lineari

Ø Attualmente più diffusi e utilizzati poiché megliosoddisfano le proprietà di un “buon” generatore,si basano su una relazione congruenziale lineareovvero sull'operazione di modulo di un intero m

Esempio: 13mod13 =

43

13= col resto di 1, cioè: 14313 +×=

y = x mod m

Ci si attende che il resto di una divisione sia sufficientemente irregolare da simulare la casualità



q Metodo congruenziale lineare

( ) mcaxx ii mod1 +=+

se è l’i-esimo valore generato, allora il successivo valore sarà:ix

dove:a: moltiplicatore

c: incremento

m: modulo

sono interi non negativi

Come seme (il valore iniziale ) viene utilizzato:

Ø l’ultimo n.ro generato da una sequenza precedente

Ø data/ora corrente

0x

mx

u ii

11

++ =

è il numero pseudocasuale generato



Esempio: 30 === xca 5=m

( )1 3 3 3 mod5 2x = ⋅ + = 4.052

1 ==u

2

40.8

5u = =( )2 3 2 3 mod5 4x = ⋅ + =

30

05

u = =( )3 3 4 3 mod5 0x = ⋅ + =

43

0.65

u = =( )4 3 0 3 mod5 3x = ⋅ + =

( ) mcaxx ii mod1 +=+ mx

u ii

11

++ =



Lo svantaggio della periodicità è ineliminabile ma può essere controllato mediante opportuna scelta delle costanti (m, a, c).Vi sono teoremi che garantiscono periodicità sufficientemente lunga

Ø m è il più grande numero intero che la macchina consente sì da garantire la migliore approssimazione dei reali. In genere è scelto primo o potenza di 2

Ø c è coprimo con m: qualsiasi se m primo, dispari se Ø a è della forma con

2m β=12 +r 2≥r

31 162 3 , 2 3 , 0m a c= + = + =Le macchine IBM 360/370 avevano un “cattivo” generatore, con:

Le macchine attuali hanno “buoni” generatori con:35 72 , 2 1, 1m a c= = + =


Generazione di numeri p.c. da v.a. qualsiasi

Finora abbiamo visto come sia possibile generare sequenze di numeri pseudocasuali provenienti da v.a. con distribuzione uniforme continua nell'intervallo (0,1).

Adesso vedremo alcuni metodi che consentono di ricavare sequenze provenienti esattamente o approssimativamente da leggi qualsiasi a partire sempre dalla distribuzione U(0,1).

In sostanza, tramite questi metodi, saremo in grado di campionare da una qualsiasi distribuzione avendo disponibile sulnostro calcolatore semplicemente il generatore uniforme.


Trasformazione integrale di probabilità

Se X è una v.a. continua la cui funzione di ripartizione è

strettamente crescente, allora la v.a. , trasformata

della X tramite la sua f.r., ha distribuzione uniforme in (0,1).

Equivalentemente

Se Y ha distribuzione uniforme in (0,1) allora la v.a.

ha proprio come funzione di ripartizione.

La trasformazione si chiama:

trasformazione integrale di probabilità

( )XF i( )XY F X=

( )1X F Y−=

( )F i

( )XY F X=


Trasformazione integrale di probabilitàSe X è una v.a. continua con f.r. esplicita e invertibile allora una determinazione x della X si ottiene applicando l'inversa della sua f.r. ad una determinazione u di una Unif(0,1):

1

0

X

Y

YI

( )1X X YI F I−=

( )XY F X=

u

( )1u Xx F u−=TIP


Trasformazione integrale di probabilitàSe X è una v.a. discreta allora formalmente cade l'invertibilità. Ma se la sua f.r. è ancora esplicita, usando la definizione di quantile, possiamo ancora ottenere una determinazione x della X a partire da una determinazione u di una Unif(0,1):

( ){ } ( ){ }inf : inf :u Xx x F x u x P X x u= ≥ = ≤ ≥

1

0

X

Y

( )XY F X=

Yu I∀ ∈

1x 2 ux x=3x

( ) ( )( 1 2,Y X XI F x F x=

TIP*


Uniforme continua in (a,b)

( )b,aUX ≈Sia ( )X

x au F x

b a−

= =−

con f.r. esplicita e invertibile

( ) ( )1Xx F u u b a a−= = − +Allora costruiamo l'inversa:

( )0,1x X u U∈ ⇔ ∈Grazie al TIP:

tramite la quale, dato un

numero pseudocasuale u,

ricaviamo immediatamente

un numero pseudocasuale x 1

0

X

U

a b

( )ux u b a a= − +


Sia X una v.c. continua con f.d. allora:( ) ≤≤

=ϕaltrove

xxxX 0

102

( )2

2

00

22

2

xx

X

tF x t dt x= = =∫

Otteniamo così la relazione:2u x=

che invertita consente di ottenere una determinazione pseudocasuale x di X:

Distribuzione continua generica

ux u=


Allora:

0 1 , 1p q p≤ ≤ = −

( ) [ ) ( ) [ ) ( )0,1 1,XF x q I x I x+∞= ⋅ +

L’algoritmo di generazione èimmediato ed esatto :

Sia u un numero pseudocasuale

Bernoulliana

Sia ( )X B p≈ con

q

1

allora

0 01 1u

se u qx

se q u≤ <

= ≤ <0

1

p


Sia ( ),X Bin n p≈

( )0

10

0 00 11 2

1

xk n k

Xk

xx

nF x p q xx

k

x n

αα

−

=

< ≤ < = = ∀ ∈≤ <

≥

∑ ¡L L

Allora:

Supponendo di aver generato un numero u tale che:

( ){ }inf : 1u Xx x F x u= ≥ =

Binomiale con TIP

0 1uα α< <

ricaveremo la nostra determinazione pseudocasuale x di X con la relazione:


Un metodo meno costoso ma ancora esatto consiste nel "simulare" n esperimenti di Bernoulli e contare gli esiti favorevoli sfruttando la proprietà:

Dunque, date n determinazioni indipendenti , 1,...,ix i n=

( )iX B p≈

1

n

u ii

x x=

= ∑

( )1

,n

ii

X X Bin n p=

= ≈∑

la loro somma X si distribuisce come una binomiale:

Binomiale per simulazione

( ) 1, ,iX B p i n≈ = …Date n bernoulliane i.i.d.

provenienti da altrettante v.a. bernoulliane ricaveremo la nostra determinazione pseudocasuale x di X con la relazione:


I metodi appena visti possono rivelarsi costosi per n grande.In questo caso conviene ricorrere all'approssimazione normale garantita dal Teorema Limite Centrale.

( )( )0.50,1

1

d

n n

X npZ N

np p →∞

− += →

−

Correzione di continuità( ),X Bin n p≈

Dunque, se z è una determinazione da una Normale Standard, otteniamo una determinazione approssimata da X con la relazione:

( )( )max 0, 0.5 1zx np z np p = − + + −

L’approssimazione Normale è accettabile per: 10 10np nq> ∧ >

Binomiale con approssimazione Normale


( )X G p≈Sia con

Allora la sua f.r. sarà: ( ) 1 xXF x q x= − ∈¥

Invertendo la relazione:ln1 x qu e= −

ricaviamo: ( )ln 1lnu

ux

q−

=

Geometrica

in ¥

0 1 , 1p q p≤ ≤ = −


( )X Exp λ≈Sia con Allora la sua f.r. sarà:

( ) ( ) [ ) ( )0,1 xXF x e I x−

+∞= −

0λ >

Invertendo la relazione: 1 xu e−= −

ricaviamo: ( )ln 1u

ux

λ−

= −

Esponenziale

in+¡


( )X P λ≈ 0>λcon( ) ,...1,0,!

=λ

= λ− xex

xfx

X

E’ semplice generare da X con TIP* (vedi binomiale) o, quandoè grande (>10), tramite approssimazione normale:

Poissoniana

λ

( )1,05.0

NX

Zd

∞→λλ →

λ+λ−

=

da cui, se z è una determinazione dalla N(0, 1) allora:

( )max 0, 0.5zx zλ λ = − + +

è un’approssimazione di una determinazione da ( )X P λ≈


Analogamente a quanto visto per la binomiale, anche nel caso della poissoniana per valori piccoli del parametro, conviene generare i tempi di arrivo, sommarli tra loro fino a che non eccedono l'unità e infine contare quanti ne sono stati generati.A garanzia del metodo interviene la proprietà:

In un processo di Poisson su un intervallo temporale unitario, iltempo di attesa tra due eventi si distribuisce come una v.a. Esponenziale di ugual parametro.

Il metodo (la cui dimostrazione omettiamo) richiede di:

Poissoniana per simulazione

1. definire k=0 e s=12. generare un numero pseudocasuale3. porre4. se altrimenti si incrementa k e si torna a 2.us e x kλ−< ⇒ =

kuks s u= ⋅


Osservando che se Z è una normale standard:

Normale non standard

( )0,1Z N∼allora: ( )2,X Z Nσ µ µ σ= + ∼desumiamo subito che, nota una determinazione pseudocasuale z proveniente da una normale standard, la corrispondente determinazione pseudocasuale proveniente da unasi ricava dalla relazione:

( )2,X N µ σ∼

zx zσ µ= +


Ø Metodo basato sul Teorema Centrale Limite

Data la successione di v.c. i.i.d. con { }i iX

∈¥

iEX µ= 2iVarX σ=e vale: ( )1,01 N

n

nX

Zn

d

n

ii

n∞→

= →σ

µ−

=∑

Per n finito ma sufficientemente elevato la v.a. nZè asintoticamente una ( )1,0N

Siano allora: indipendenti. Avremo che: ( )0,1 1, ,iU U i n=∼ …( )

21/2 , 1/12 , 0,1

12

dii

i i n n

U nEU VarU Z N

n →∞

−= = = →∑

Ø particolareØ approssimato

Normale standard tramite TLCPoiché la f.r. di una non è data in forma esplicita nonè possibile ricorrere al metodo della TIP

( )1,0NZ ≈


Dunque, dati 12 numeri pseudocasuali 12,...,1, =iui

Così, per n=12, si ha:12

121

6ii

Z U=

= −∑

è una (buona) approssimazione di una determinazione proveniente da una normale standard

12

1

6u ii

z u=

= −∑

Osservazioni:• Metodo piuttosto lento (servono ben 12 pseudocasuali)• Z è troncata fuori da (-6, +6)• Il metodo fornisce soluzioni approssimate ma accettabili

Normale standard tramite TLC


Ø particolareØ esatto

Metodo di G.E.P. Box e M.E. Mullero delle coordinate polari - 1958

Normale standard tramite Box-Muller

Consideriamo la gaussiana standard bivariata nelle coordinate polari :( ),r θ

( ) ( ) ( )2 / 2,2

rrg r drd e drd p r q drdθ θ θ θ θ

π−= =

e calcoliamo le f.r. marginali del modulo e dell'angolo del vettore polare:

( ) ( )

( ) ( )

2 / 21

0

20

1

2

rru P r p t dt e

u Q q t dtθ θ

θπ

−= = = −

= = =

∫

∫


Normale standard tramite Box-Muller

Poiché entrambe le f.r. sono analiticamente invertibili, possiamo applicare la TIP e ricavare facilmente:

1

2

2 ln

2

r u

uθ π

= −

=

Ora torniamo alle coordinate cartesiane ricavando una coppia di determinazioni pseudocasuali normali standard indipendenti a partire da una coppia di determinazioni pseudocasuali uniformi:

( )( )

1 1 2

2 1 2

cos 2ln cos 2

sin 2ln sin 2

z r u u

z r u u

θ π

θ π

= = −

= = −


Metodo Box-Muller ottimizzatoPoiché il computo delle funzioni trigonometriche è estremamente dispendioso si può ricorrere al seguente ingegnoso stratagemma.

Se generiamo casualmente un punto (x,y) all'interno del cerchio goniometrico è possibile dimostrare che:

( )

( )

2 21 2 2

2 22 2 2

2 ln

2 ln

xz x y

x y

yz x y

x y

= − + + = − + +

( ) ( ) ( )2 2

2 2 2 20,1 , cos 2 0,1 , sin 2 0,1

X YX Y U U U

X Y X Yπ π+

+ +∼ ∼ ∼

Allora la precedente relazione può essere riscritta in termini di x e y senza l'uso delle funzioni trigonometriche:

Date post:	20-Feb-2019
Category:	Documents
Upload:	dodiep
View:	216 times
Download:	0 times

Numeri casuali -...

Documents