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