1 Importance sampling: metodo analitico.
Di seguito proviamo a verificare graficamente le formule trovate analiticamente nelle dispense per ottenere una sequenza di numeri casuali che segua una certa distribuzione . Il primo esempio è quello della seguente funzione peso , che integrata e invertendo ci fornisce la relazione , con y numeri casuali distribuiti uniformemente mentre le saranno invece distribute come la funzione peso . Mediante le poche righe di codice seguente proviamo a verificare che sia effettivamente cosi:
1 2 3 4 5 6 7 8 9 10 11 | N=10000
y=rand(N)
x=2-sqrt(4-3*y)
w=(4-2*x)/3
hist(x,10,normed=True)
plot(x,(4-2*x)/3,'s')
#sorting
x.sort() #sul posto, x viene sostituito da quello ordinato
x_sorted = sort(x) #il risultato va memorizzato su un altro array
|
Come notiamo dal grafico, la distribuzione dei valori dell'array x e la funzione peso analitica sono in buon accordo. L'istogramma è stato realizzato mediante la funzione hist(x,10,normed=True): il primo argomento è l'array d'interesse, il secondo è il numero di bin, intervalli, all'interno dei quali contare quanti valori di x si hanno, con il terzo chiediamo che l'istogramma sia normalizzato. Il plot della funzione lo si fa per punti (notare il 's') perchè la sequenza x non è ordinata. L'ordinamento è utile in qualche caso e lo si fa mediante la funzione sort().
Ora che abbiamo generato la nostra sequenza random distribuita come , possiamo calcolare la somma 182 e confrontarla con la somma fatta usando una distribuzione uniforme, in questo modo:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | fx=1/(1+x**2)
f=1/(1+y**2)
I_esatto = pi/4
I_imp = sum(fx/w)/N
I_unif = sum(f)/N
err_imp=log10(abs(I_imp - I_esatto))
err_unif=log10(abs(I_unif - I_esatto))
print 'Importance sampling: ' + str(I_imp) + ' Log10 error: ' + str(err_imp)
print 'Uniform sampling: ' + str(I_unif) + ' Log10 error: ' + str(err_unif)
#Output
Importance sampling: 0.785375419649 Log error: -4.64313796351
Uniform sampling: 0.787076269334 Log error: -2.77518062628
|
Abbiamo creato un array fx che contiente i valori della funzione campionata secondo i valori di x, da dividere poi per la funzione w. L'array f invece contiene i valori della stessa funzione campionata in maniera uniforme pero usando la sequenza random y. Come si può vedere dall'output, l'integrale calcolato usando l'importance sampling è piu vicino a quello vero, avendo un errore di due ordini di grandezza inferiore rispetto all'altro caso.
Volendo riprodurre il grafico dell'andamento dell'errore in funzione del numero di punti random di campionamento, mostrato nelle dispense in fig. 46, è necessario eseguire le righe di codice viste all'interno di un ciclo for nel quale modifichiamo il valore della variabile N. Questo però non basta a causa delle normali fluttuazioni nei valori dell'integrale, percui non otterremo degli andamenti quali lineari come in figura ma delle nuvole di punti. Per ridurre le fluttuazioni, dobbiamo eseguire piu volte il calcolo dell'integrale e fare la media per ogni valore di N; ci servirà perciò un ulteriore ciclo for. Lasciamo come esercizio la scrittura di questo codice e la produzione del grafico in figura 46. Inoltre, lasciamo a voi la verifica delle altre distribuzioni ottenute per via analitica, come la esponenziale (formula 185), uniforme su intervallo generico (186) e gaussiana (191).
2 Importance sampling: algoritmo Metropolis.
Quando non si conosce la forma analitica della funzione distribuzione che si vuole genenerare o non si può applicare il metodo analitico, si ricorre all'algoritmo Metropolis, uno dei più importanti algoritmi mai inventati.
- si sceglie una posizione di prova , dove è la posizione precedente e è un numero random nell'intervallo ;
- si calcola la probabilità ;
- se , si accetta la nuova posizione, quindi ;
- se , si genera un numero random tra 0 e 1;
Questi semplici passaggi, che si traducono in poche righe di codice, ci permettono di avere una sequenza di numeri random distribuiti come la funzione . Vediamo di seguito una implementazione che segue i passi appena visti:
1 2 3 4 5 6 7 8 9 | xt=xi+delta*(2*rand()-1)
w=p(xt)/p(xi)
if w >=1:
xi=xt
elif w < 1:
r=rand()
if r < w:
xi=xt
|
Ovviamente, per poter funzionare vanno definite le variabili delta, xi e la funzione p(x). Inoltre queste righe genererebbero un solo numero casuale, xt appunto, a partire dal precedente xi. Prima di completare la nostra routine, vediamo di seguito una versione più efficente:
1 2 3 4 5 | xt=xi+delta*(2*rand()-1)
w=p(xt)/p(xi)
r=rand()
if w > r:
xi=xt
|
Può sembrare strano ma queste righe fanno esattamente la stessa cosa di quelle scritte in precedenza. Il trucco sta nel fatto che il numero random è sempre compreso tra 0 e 1, perciò quando sarà sempre verificata anche la condizione , quindi la mossa viene sempre accettata; l'altra condizione è che , che è quella che usiamo e che di fatto racchiude in se entrambe i due casi.
Quello che dovremmo aggiungere è un ciclo sul numero di mosse accettate e una lista/array che conservi tutte le mosse accettate, che costituiscono la sequenza di numeri random che stiamo cercando di generare. Vanno inoltre definite alcune parametri come delta e la posizione iniziale e la funzione peso che vogliamo usare. Completiamo quindi il nostro codice come segue.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | def p(x):
return 2*exp(-(x)**2)+exp(-(x+3)**2)+exp(-(x-3)**2)
delta=3
x=[0.]
N=5000
naccept=0
counter=0
while naccept < N:
counter+=1
xt=x[-1]+delta*(2*rand()-1)
w=p(xt)/p(x[-1])
r=rand()
if w > r:
x.append(xt)
naccept+=1
print naccept/counter
|
Per poter salvare tutte le posizioni accettate usiamo una lista x di cui poniamo a 0 il primo elemento che rappresenta la posizione di partenza. All'interno del ciclo while usiamo x[-1] per leggere l'ultima posizione registrata/accettata e usiamo la funzione append() per aggiungere alla lista la posizione xt in caso questa sia stata accettata. Il ciclo è eseguito sintanto che la variabile naccept raggiunge il valore desiderato, 1000 in questo caso, che sarà anche la lunghezza finale della lista x. Per trovare un valore accettabile di delta è utile controllare il rapporto tra tutti i passi generati e quelli accettati, come facciamo nell'ultima riga: questo rapporto normalmente dovrebbe essere minore di 0.5 per avere delle sequenze accettabili. Ora vediamo se la nostra sequenza random x è distribuita come vogliamo, cioè come la funzione.
1 | hist(x,100,normed=True)
|
Cosi creiamo direttametne un istogramma normalizzato con 100 intervalli. Altrimenti, per fare un grafico per punti:
1 2 | h,b=histogram(x,100,normed=True)
plot(b[:-1],h)
|
Come si nota, anche con una sola esecuzione del programma e usando non tantissimi passi, l'istogramma è piuttosto simile alla funzione analitica che volevamo riprodurre. Se volessimo avere una maggiore precisione dovremmo aumentare la lunghezza della sequenza e fare delle medie su più istogrammi.
Quando si a che fare con funzioni periodiche, è bene limitare la generazione di nuove posizioni all'interno di un intervallo stabilito. Altrimenti il walker andrà in giro per tutto lo spazio cercando di riprodurre tutte le copie periodiche della funzione. Per esempio volendo riprodurre la funzione di hanning con periodo 10, senza limiti e un delta di 1 e 5, otteniamo:
Nella prima figura si riconoscono 2 periodi positivi della funzione, non riprodotti al meglio visto che hanno diversa altezza. Molte più copie si hanno se il delta è 5, ma l'accuratezza diminuisce ulteriormente. E' perciò doveroso porre un limite alle posizioni generate per replicare accuratamente l'intervallo che si desidera. Potremo modificare il ciclo while, introducendo una condizione if sui valori della xt, come segue:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | def p(x):
return 0.5-0.5*cos(2*pi*x/10)
delta=1.
N=50000
x=[5.]
naccept=0
counter=0.
while naccept < N:
xt=x[-1]+delta*(2*rand()-1)
counter+=1
if xt <10 and xt > 0:
w=p(xt)/p(x[-1])
r=rand()
if w > r:
x.append(xt)
naccept+=1
print naccept/counter
hist(x,100,normed=True)
xx=arange(0,10,0.1)
plot(xx,p(xx)/4.5)
|
In questo modo escludiamo le posizioni che sono fuori dall'intervallo , in questo modo tutti le N posizioni generate saranno al suo interno. Facciamo notare come la posizione iniziale sia in questo caso 5, sempre in corrispondenza del massimo della funzione. Con le ultime due righe facciamo un grafico anche della funzione analitica (riscalata). Questo il risultato.