Importance sampling: metodo analitico.

 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

 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

Importance sampling: algoritmo Metropolis.

  1. si sceglie una posizione di prova x t = x i + δ i , dove x i è la posizione precedente e δ i è un numero random nell'intervallo [ - δ , δ ] ;
  2. si calcola la probabilità w=p( x t )/p( x i ) ;
  3. se w1 , si accetta la nuova posizione, quindi x i+1 = x t ;
  4. se w<1 , si genera un numero random tra 0 e 1;
    1. se r<w , si accetta la nuova posizione, quindi x i+1 = x t ;
    2. altrimenti la posizione viene rifiutata
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 &lt; 1:
	r=rand()
	if r &lt; w:
		xi=xt

1
2
3
4
5
xt=xi+delta*(2*rand()-1)
w=p(xt)/p(xi)
r=rand()
if w > r:
	xi=xt

 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 &lt; 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

1
hist(x,100,normed=True)

1
2
h,b=histogram(x,100,normed=True)
plot(b[:-1],h)

 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 &lt; N:
	xt=x[-1]+delta*(2*rand()-1)
	counter+=1
	if xt &lt;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)

Simulated annealing