Metodi Monte Carlo: numeri pseudo casuali e importance sampling

Numeri Pseudo Casuali

Integrazione Monte Carlo Numeri pseudo casuali Prima di introdurre i vari metodi, chiariamo cosa funziona la simulazione di variabili aleatorie. Il punto fondamentale è che un computer non può simulare qualcosa che sia "random", ma solo valori pseudo-random, quindi deterministici (P.S .: esistono siti che vendono numeri "veramente" casuali, per esempio https: //www.random.org) Partiamo dal caso uniforme, possiamo dire che Definizione - uniform Pseudo Random generator Un generatore di variabili pseudo-casuali U(0, 1), è un algoritmo che partendo da un punto iniziale u0 (seed) e una trasformazione D(), produce una sequenza di valori (u1, ... , Un), con ui = D(ui-1) € [0, 1]. Il vettore (u1, ... , Un) riproduce le "caratteristiche" di un vero vettore di realizzazioni di variabili aleatoria 1,. , Vn) iid da U(0,1).

Integrazione Monte Carlo

Osservazioni sui numeri pseudo casuali

• L'algoritmo è deterministico, e conoscendo il seme (seed) e la funzione D(), possiamo ricostruire il campione pseudo-casuale;
• Con "riproduce le caratteristiche" si intende che nessun test del tipo Ho : 01, . . , Un id U(0,1) venga rifiutato.
• Le sequenze sono periodiche.

L'idea è quindi che la sequenza di numeri pseudo casuali sono indistinguibili da numeri casuali. In R possiamo fare quest'esempio. set . seed (100) runif (10,0,1) set. seed (100) runif (10,0,1) runif (10,0,1)

Integrazione Monte Carlo: Generatore U(0,1)

Nel codice sopra, i primi 2 runif daranno lo stesso risultato (perchè risettate il seed), nel terzo i risultati sono diversi. Un semplice esempio di generatore pseudo casuale U (0, 1) si può ottenere creando le variabili di = (adi-1 + c) mod m dove 0 ≤ c ≤ m (incremento) e 0 < a < m (moltiplicatore), e dopo prendendo Ui = di / m. I valori m, a e c, vanno scelti in modo da non creare loop.

Integrazione Monte Carlo: Scopo

Integrazione In molti studi Monte Carlo, lo scopo è studiare delle caratteristiche delle variabili aleatorie X E X' nella forma E(h(X) =/ h(x)f(x)dx se i dati sono continui, oppure E(h(X)) => xEx h(x)P(X =x) se discreti e E(h(X)) = > XEXD h(x)P(X=x) + / Xc h(x)f(x)dx se la variabile è mista con X'D parte del dominio discreto e Ac continuo con XD UXC = X. Oppure, per scriverli in una forma unificata E(h(X)) = |h(x)f(x)d)(x) X con dx(x)

Integrazione Monte Carlo: Misure

• misura di Lebesgue: d)(x) = dx nel caso continuo
• misura di conteggio nel caso discreto.
• misura ibrida in caso di variabili miste

Esempi di h(X)

• h(X) = X - Media;
• h(X) = X2 momento secondo . ..

Integrazione Monte Carlo: Convergenza

Dato un campione di variabili identicamente distribuite, ma non necessariamente indipendenti (hanno quindi la stessa marginale) (X1, ... , Xn), sappiamo che hn = Xi=1h(Xi) n converge quasi certamente a E(h(X)) per la legge dei grandi numeri => Questo è vero anche per h() diversa dall'identità' visto che se definiamo Y = h(x), e indichiamo con fy (y) la sua densità, allora Ix h(x)f(x)dx(x)= |yfy(y)d)(y) ~ Y 2 Se h2(X) ha attesa finita, possiamo calcolarne anche la varianza Un = var (hn) = var n Zi 1h(Xi) i=1 n ◆ = 1 n2 n i=1 var (h(Xi)) + i=1 j=1 >> cov (h(Xi), h(Xj)) n ))!

Integrazione Monte Carlo: Varianza

Nel caso in cui le variabili siano indipendenti abbiamo che var(hn) = var (h(Xi)) n = 1 n Jx (h(x)- E(h(X)))2f(x)dx(x) che può essere approssimato da Un = p2 2(h(xi) - hn)2 n.2 n i=1 Attenzione: i) Lo stimatore della varianza non è uno stimatore Monte Carlo. ii) Si può dimostrare che è distorto iii) Una versione non distorta dello stimatore della varianza si trova sostituendo(n - 1) a n nel denominatore Naturalmente più piccola è var(hn) e migliore sarà l'approssimazione di E(h(X)). Da ora in poi, a meno che non sia detto esplicitamente, assumeremo che le variabili siano anche indipendenti

Integrazione Monte Carlo: Funzioni Multivariabili

I Monte Carlo vale anche per funzioni di più variabile: E(h(X,Y) = h(x,y)f(x,y)d)(y)d)(x) ~ n n i=1 h(xi, yi) e anche nel caso in cui h(x, y) = h*(x): E(h*(X)) = | h*(x)f(x)d)(x) = X n i=1 h (xi, yi) n = >71 h* (xi) i=1 h(x,y)f(x,y)d)(y)d)(x) ~ XIV n

Integrazione Monte Carlo: Funzione di Ripartizione

Funzione di ripartizione La funzione di ripartizione per una variabile può essere scritta come F(a) = [° O f(x)dx(x) che non è nella forma che ci permette di usare un metodo Monte Carlo. Possiamo però introdurre la funzione 1a (x) che assume valore 1 se x < a, altrimenti zero, e usare f(x)dx(x) = -00 ra 1 -00 00 1a(x)f(x)d)(x)=E(1a(X)) ~ n Xi=1 1a (Xi) i=1 dove naturalmente Xi ~ F e in questo caso h(X) = 1a (X). Possiamo quindi stimare l'intera funzione di ripartizione con Ê(a) = i=1 2,1 1a (Xi) n ~ E(1a(X) = / 00 1a (x)f(x)dx = 1 a -00 f(x)dx = F(a) Ê(a) è una stima della funzione di ripartizione

Integrazione Monte Carlo: ECDF

1.0 0.8 0.6 ecdf 0.4 True n=10 0.2 n=50 n=100 0.0 0.0 0.2 0.4 0.6 0.8 1.0 X Figure: ecdf di una B(1,3) per diversi valori di n

Integrazione Monte Carlo: Codice ECDF

Code: Codice della figura # simuliamo dalla beta x1 = rbeta (10,1,3) x2 = rbeta (50,1,3) x3 = rbeta (100,1,3) # cacoliamo le funzioni di ripartizioni empiriche su xseq xseq= seq(0,1, by=0.001) ec1 = ecdf (x1) (xseq) ec2 = ecdf (x2) (xseq) ec3 = ecdf (x3) (xseq) plot (xseq, ec1, type="s", lwd=2, col=2, xlab="X", ylab="ecdf") lines (xseq,ec2, type="s", col=3, lwd=2) lines (xseq,ec3, type="s", col=4, lwd=2) lines (xseq, pbeta(xseq,1,3), type="s", col=1, lwd=2) legend ("bottomright ", c ( "True " , "n=10" , "n=50" , "n=100") , col=1:4, lwd=3, cex=2 )

Integrazione Monte Carlo: Quantili Empirici

Come si vede dalla figura la funzione è a scalini, con scalini di altezza 1/n e per Xi-1 < a ≤ Xi è costante. Dalla funzione di ripartizione stimata (chiamata anche funzione di ripartizione empirica), si possono calcolare i quantili empirici. Per esempio il quantile empirico di livello p è min{x; : i = 1, ... , n and F(x) = p}

Integrazione Monte Carlo: Stima di Densità

Stima di densità Un caso interessante è quando non sappiamo come è fatta f (x) (forma funzionale), però conosciamo f (x)y), dove Y ~ G è definita su ), e sappiamo campionare da G. Possiamo calcolare la densità f (x) in un punto specifico x *: f(x*) = f(x*|y)g(y)dx(y) Y e se indichiamo con Y1, ... , Yn campioni iid da G, e assumiamo h(X) = f(x*|y), allora _,_1f(x*|Yi) i=1 ⇡ Z Y f(x⇤|y)g(y)d n

Integrazione Monte Carlo: Densità Marginale

Esempio di densità marginale

Ipotizziamo che Y ~ Exp(X), e X|Y = y ~ G(expy, 1), come è fatta la densità marginale di X ?. Questo sopra è un semplice caso in cui non siamo in grado di dire qualcosa su una distribuzione e possiamo usare l'approssimazione MC 0.4 - d 0.3 fx_stima 0.2 0.1 0.0 7 - - O 0 2 4 6 8 10 a Figure: Esempio MC

Integrazione Monte Carlo: Codice Densità Marginale

Code: Codice della figura n = 1000 # valori di x su cui calcolare la densita' a = seq(0,10,by=0.1) lambda = 2 fx_stima = c () y = rexp (n, lambda) for(i in 1: length (a)) { # densita ' condizionata fzgiveny = dgamma(a[i],exp(y),1) # stima di fx nel punto a[i] fx_stima [i] = sum(fzgiveny) /n } plot (a, fx_stima) lines(a,fx_stima, col=2) Naturalmente non sappiamo determinare la forma funzionale, ma possiamo trovarla utilizzando qualche metodo interpolante, o usando funzioni a tratti.

Rao-Blackwellization

Prendiamo la seguente relazione E(h(X) = / h(x)fx(x)d)(x) = h(x)f(x|y)f(y)d)(x) dx(y) = Y Z Y Z X X h(x)f(x|y)d)(x) )f(y)d)(y) e quindi Ex (h(X)) = Y Ex(h(X)|y)f(y)d)(y) quindi questi valori attesi potremmo calcolarli usando i campioni di X, con Zi=1 h(Xi) i=1 n ~ E(h(X)) o con quelli di Y, con 27-1 Ex (h(X)|yi) i=1 ⇡ Ex(h(X)|y)f(y)d)(y) n

Rao-Blackwellization: Varianza

se siamo in grado di calcolare E(h(X)|Yi). Questa eguaglianza tra i valori attesi è la classica legge E(X) = Ey(Ex(X|Y)) Assumendo indipendenza, abbiamo che Zi=1 var x (h(Xi)) varx = i=1 n2 = e vary In i=1 Ex (h(X)|Yi) ◆ = Zi 1 vary (Ex (h(X)|Y;)) n2 Sebbene i due valori attesi stimano/calcolano la stessa cosa, abbiamo che varx (h(X)) = vary (Ex (h(X)|Y)) + Ey (varx (h(X)|Y))) ≥ vary(Ex(h(X)|Y)) (legge della varianza totale) Abbiamo quindi che varx i=1 Zi=1h(Xi) n ◆ z vary i=1 Ex (h(X)|Yi) n ) n Xi-1 h(Xi) n ◆

Rao-Blackwellization: Stimatore

e quindi l'approssimazione Ei=1 n n Ex (h(X)|Yi) è migliore in termine di varianza. Lo stimatore che utilizza la distribuzione condizionata per calcolare il valore atteso marginale di h(X), si chiama Rao-blackwell estimator.

Importance Sampling

Importance sampling Ritorniamo al problema del calcolo di attese E(h(x) = / h(x)f(x)dx(x) Questo integrale può essere scritto come E(h(X) = h(x)f(x)dx(x) = g(x) h(x)f(x) g(x)dx(x) dove g(x) è una qualsiasi funzione per cui g(x) > 0 per ogni x € X' tale che f (x) > 0. Nel caso in cui g(x) è una densità abbiamo che Z x h(x)f(x) g(x) g(x)dx(x) = Eg h(X)f(X) g(X) ◆ Come per l'esempio precedente, X h(x)f(x)dx(x)

Importance Sampling: Stimatori

e X h(x)f(x) g(x) g(x)dx(x) calcolano la stessa cosa, Eg h(X)f(X) g(X) ◆ = Ef(h(X)) Abbiamo quindi due possibili stimatori . con Xi ~ f Ei=1h(xi) i=1 ~ Ef(h(X)) n · con Xi ~ g Li=1 g(xi) in h(xi)f (xi) n ~ Eg h(X)f(X) g(X) ◆

Importance Sampling: Varianza

Possiamo vedere quale stimatore ha varianza minore, se assumiamo indipendenza, e calcoliamo varf(h(X) = E(h(X)2) - E2(h(X) = h2(x)f(x)d)(x)-E?(h(X)) e varg g(X) h(x)f(x) = Jx g2 (x) h2(x)f2(x) g(x)d)(x)- E}(h(X)) dove, nella seconda varianza, abbiamo usato la relazione Eq 12 h(X)f(X) g(X) ◆ = E3(h(X)) Abbiamo allora che l'importance sampling è migliore se varf (h(X)) - varg h(X)f(X) g(X) ) >0 Abbiamo che var f (h(X)) - varg g(X) h(x)f(x) = [ (1 - f(x) g(x) ◆ h2(x)f(x)dx(x)

Importance Sampling: Scelta di g(x)

visto che f (x) e g(x) sono densità, il rapporto f(x)/g(x) non può sempre essere <1 o >1. Quindi, per avere una differenza di varianze positiva (i.e., importance sampling migliore del vanilla), dobbiamo avere che · f(x)/g(x) > 1 se h2(x)f(x) è piccolo; 1 - f(x) g(x) negativo; · f(x)/g(x) < 1 se h2(x)f(x) è grande; 1 - g(x) f(x) positivo. quindi g(x) deve essere elevata (i.e., importante), per i punti dove h2 (x)f(x) è elevata. Si può dimostrare che la scelta migliore per g(x) è g(x) = |h(x)|f(x) Jx |h(x)|f(x)d)(x)

Simulazioni

Simulazione Se sappiamo generare da una distribuzione, con trasformazioni di variabili possiamo ottenere campioni da altre distribuzioni.

Esempio di simulazione

Se U ~U(0,1), allora X =- logU~Exp(1) e ~~ Exp(A) Soluzione: possiamo usare la regola della trasformazione di variabili: fx(x) = fu(u) du dx =- exp(-x)| = exp(-x) e fY(y) = fx (x) dy dx = exp(-Xy)|X|=Xexp(-Xy)