simulacion_montecarlo_pi

simulacion_montecarlo_pi

#!/usr/bin/env python
#-*- coding:utf-8 -*-

from __future__ import division
import numpy as np
from pylab import *

pi_hat = []
nn = 1000 # numero de veces que repetimos la simulación

for i in np.arange(nn):

	n = 1000 # número de valores aleatorios
	u = np.random.uniform(0,1,n) # valores aleatorios de un eje
	v = np.random.uniform(0,1,n) # valores aleatorios del otro eje
	dis = np.sqrt((u-0.5) ** 2 + (v-0.5) ** 2) # distancia al centro
	dis_i = filter(lambda x: x < 0.5, dis) # filtramos los que cumplen la condición de círculo
	o = len(dis_i) / n # calculamos el área del círculo
	pi = o * 4 # aplicamos la ecuación de pi
	pi_hat += [pi] # guardamo la estimación de pi de cada simulación

print "Pi_hat = ", np.mean(pi_hat) # mostramos el valor medio de pi para las nn estimaciones 
print "Pi real = 3.14159265359" # el valor real de pi
print "ECM =", np.mean(np.subtract(3.14159265359,pi_hat) ** 2)

# Funciones auxiliares para el grafico
# necesitamos solo los valores aleatorios que cumplen la condicion de circunferencia
u_f = [] 
v_f = []

for i,e in enumerate(u):
    if np.sqrt((u[i]-0.5 )** 2 + (v[i]-0.5) ** 2) < 0.5:
        u_f += [u[i]]
        v_f += [v[i]]

plt.figure(figsize=(6,12))

ax1 = subplot(211) # gráfico de la simulacion de montecarlo
ax1.axhline(0.5, color="grey", alpha=0.5)
ax1.axvline(0.5, color="grey", alpha=0.5)
text(0.2, 1.25, u"Simulación Montecarlo", fontsize=14)
plt.scatter(u,v)
plt.scatter(u_f,v_f, linewidth=0.08, color="green")

ax2 = subplot(212) # gráfico de las estimaciones de pi
plt.scatter(range(len(pi_hat)),pi_hat)
ax2.axhline(y=3.14159265359, color="grey", alpha=0.5, linewidth=2)
text(len(pi_hat)*1.01, 3.14159265359, "Pi = 3.1415...")
ax2.set_title(u"Estimación de Pi", fontsize=14)

#plt.savefig("simulacion_pi_montercalo.png", dpi=200) # quitar el \# para guardar la imagen.