-
Notifications
You must be signed in to change notification settings - Fork 1
/
3.a. metropolis.py
68 lines (47 loc) · 1.19 KB
/
3.a. metropolis.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
Este script genera una muestra aleatoria con
PDF \propto x / (1 + x**2)**
usando el algoritmo de Metropolis
'''
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import uniform
from numpy.random import seed
def W(x):
return x / (1 + x**2)**2
# def W(x):
# return 3.5 * np.exp(-(x-3.)**2/3.) + 2 * np.exp(-(x+1.5)**2/0.5)
x = np.linspace(0, 8, 10000)
plt.figure(1)
plt.clf()
# plt.plot(x, W(x) / 13.251318549)
plt.plot(x, W(x) / 0.49)
plt.xlim(-2, 10)
plt.xlabel(r"$\{x_i\}$", fontsize=18)
plt.ylabel(r"$\propto W(x)$", fontsize=18)
# plt.show()
# plt.draw()
from scipy.integrate import trapz
print "La integral de W(x) es = {}".format(trapz(W(x), x))
# Metropolis
seed(123)
d = 1.
N_muestra = int(1e5)
x0 = 8
# un paso
def paso_metropolis(x0):
r = uniform(low=-1, high=1)
xp = x0 + d * r
if W(xp) / W(x0) > uniform(0, 1):
x0 = xp
return x0
muestra = np.zeros(N_muestra)
muestra[0] = x0
for i in range(1, N_muestra):
muestra[i] = paso_metropolis(muestra[i-1])
plt.hist(muestra, bins=np.arange(-4, 8, 0.1), normed=True, histtype='step')
plt.show()
plt.draw()