import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
%matplotlib inline
Histograma y sus errores¶
El número de datos en cada bin del histograma es una variable aleatoria. La distribución de esta variable es binomial, pues cada uno de los $N$ resultados de las simulaciones tienen una cierta probabilidad $p$ de caer en el intervalo que corresponde a una clase y se cuentan la cantidad de éxitos $k$, i.e. $Binom(N, p, k)$ da la probabilidad de contar $k$ casos en $N$ casos totales. Entonces, para estimar el error en cada bin usamos la desviación estándar de la distribución binomial, $\sigma = \sqrt{N \, p \, (1-p)}$.
Si queremos representar una distribución de probabilidad, necesitamos normalizar el histograma. Para esto, dividimos tanto los bins como su error por su ancho y el numero total de datos del histograma $N$.
Sin embargo, a veces tiene sentido considerar histogramas con ancho de bin variable. Usar bins más anchos donde la densidad es menor reduce el ruido debido a la aleatoriedad del sampleo. Usar bins más angostos donde la densidad es mayor da mayor precision a la estimacion de la densidad de probabilidad.
Para estimar la probabilidad $p$ de la binomial, consideramos a los $k_i$ datos que cayeron en el bin $i$ y lo dividimos por el numero total de datos $N$. Es decir, $p_i = k_i/N$. En el caso de que los bins tengan distinto ancho $b$, hay que tenerlo en cuenta al calcular la probabilidad:
$$p_i = \frac{b_i k_i}{\sum_j b_j k_j}$$Cuando los anchos son todos iguales, se recupera la formula anterior.
def hist(data, bins=10, normed=False, **kwargs):
'''Grafica histogramas con barras de error'''
y, bins = np.histogram(data, bins=bins) #Divido los datos en bins
binwidth = (bins[1:]-bins[:-1])
N = np.sum(y) * binwidth if normed else 1
yerr = hist_err(y, binwidth=binwidth) / N #Calculo y normalizo el error
y = y / N #Normalizo los datos
f, ax = plt.subplots(1,1)
ax.bar(bins[:-1], y, width=binwidth, yerr=yerr, **kwargs)
def hist_err(data, binwidth=1):
'''Calcula el error (binomial) para un histograma'''
N = np.sum(data * binwidth)
p = data * binwidth / N
return np.sqrt(p * (1 - p) * N)
errdict = dict(ecolor='red', elinewidth=3, capsize=5, capthick=3, alpha=0.75)
formato = dict(color='r', alpha=0.5, error_kw=errdict)
Guia 3, ejercicio 4 (b)¶
Vamos a aproximar la densidad de probabilidad de Cauchy por una suma de 2 gaussianas centradas en 0 y con distintos anchos.
$$f (x) = a \, N (0, \sigma_1) + b \, N (0, \sigma_2)$$donde $a + b = 1$ porque tiene ser una probabilidad. Es decir:
$$ 1 = \int f(x) \, dx = \int a \, N (0, \sigma_1) + b \, N (0, \sigma_2) = a + b$$Si aumentamos el rango, la aproximación deja de tener validez. Esto se debe a que la gaussiana decae mucho más rápido que la Cauchy. Vemos en el segundo gráfico que para $|x|>10$ la aproximación por gaussianas ya toma valores muy pequeños.
De hecho, la varianza de la distribución de Cauchy no está definida, es infinita. Por lo que si quisieramos aproximar la densidad de Cauchy por gaussianas, estas tendrían que tener un ancho infinito.
def cauchy(x):
return st.cauchy.pdf(x)
def cauchy_apr(x):
'''Aproximacion por 2 gaussianas para Cauchy'''
return 0.5 * st.norm(scale = 0.75).pdf(x) + 0.5 * st.norm(scale = 3).pdf(x)
x1 = np.linspace(-6, 6, 100)
plt.plot(x1, cauchy(x1)/cauchy(0), label='Cauchy')
plt.plot(x1, cauchy_apr(x1)/cauchy_apr(0), label='Aprox')
plt.legend()
x2 = np.linspace(-30, 30, 1000)
plt.plot(x2, cauchy(x2)/cauchy(0), label='Cauchy')
plt.plot(x2, cauchy_apr(x2)/cauchy_apr(0), label='Aprox')
plt.ylim((0, 0.01))
plt.legend()
x1 = np.linspace(-30, 30, 1000)
plt.plot(x1, cauchy(x1)/cauchy(0), label='Cauchy')
plt.plot(x1, cauchy(x1)/cauchy(0) - cauchy_apr(x1)/cauchy_apr(0), label='Diferencia')
plt.ylim((-0.03, 0.03))
plt.legend()
Guia 3, ejercicio 9 (b)¶
Vamos a generar una variable Y con distribución exponencial a partir de una variable X con distribución uniforme.
Como la probabilidad se tiene que conservar, tenemos:
$$ f_y(y) \, dy = f_x(x) \, dx $$$$ a \, exp(-ay) \, dy = 1 \, dx $$$$ x(y) = exp(-ay) $$$$ y(x) = -\frac{log(x)}{a} $$Ancho de bins¶
Como se habló antes en la seccion de histogramas, en esta distribución puede ser conveniente usar bins de ancho variable, y aumentar el ancho cuando nos alejamos del origen, donde se espera obtener cada vez menos datos.
Segun lo comentado en clase, hay diversos criterios y ninguno es el 'correcto'. Asi que para el segundo grafico usé bins espaciados logaritmicamente despues de $y>15$.
def exp(x, a):
return a*np.exp(-a*x)
a = 0.1 #Parametro de la exponencial
x = st.uniform.rvs(size=500)
y = -np.log(x) / a
hist(y, bins=20, normed=True, **formato)
x = np.linspace(0, np.max(y), 1000)
plt.plot(x, exp(x, a), 'g', lw=3, alpha=0.75)
bins1 = np.linspace(0, 15, 5, endpoint=False)
bins2 = np.logspace(np.log(15), np.log(np.max(y)), 5, base=np.e)
bins = np.concatenate((bins1, bins2))
hist(y, bins=bins, normed=True, **formato)
x = np.linspace(0, np.max(y), 1000)
plt.plot(x, exp(x, a), 'g', lw=3, alpha=0.75)
Guia 4, ejercicio 10¶
La densidad conjunta de probabilidades es:
$$ f(r, \theta, \phi) = \frac{1}{4\pi} \, \frac{1}{1 - \pi/4} \, \frac{1}{1 + r^2} \, r^2 \, sin \, \theta \, (0 < r < r_e)$$El factor $\frac{1}{4\pi} \, \frac{1}{1 - \pi/4}$ proviene de la normalización: $\int f(r, \theta, \phi) \, dr \, d \theta \, d\phi = 1$.
Para calcular las distribuciones marginales de cada variable, integramos la probabilidad conjunta sobre el resto de las variables:
$$ f_r(r) = \int f(r, \theta, \phi) \, d \theta \, d\phi = \frac{1}{1 - \pi/4} \, \frac{r^2}{1 + r^2} $$$$ f_\theta(\theta) = \int f(r, \theta, \phi) \, dr \, d\phi = \frac{sin \, \theta}{2} $$$$ f_\phi(\phi) = \int f(r, \theta, \phi) \, dr \, d \theta = \frac{1}{2 \pi} $$Entonces, $f(r, \theta, \phi) = f_r(r) \, f_\theta(\theta) \, f_\phi(\phi)$, y por lo tanto son independientes.
Aplicando el método de Monte Carlo, generamos la distribución de las variable $r$ y $\theta$. Para la variable $\phi$ corresponde una distribución uniforme.
def fr(x):
#Distribucion de r
return 1/(1-np.pi/4) * x**2/(1+x**2)
def ft(x):
#Distribucion de theta
return np.sin(x) / 2
r = []
while len(r) < 10000:
u = st.uniform.rvs() #Uniforme entre 0 y 1
v = st.uniform.rvs() #Uniforme entre 0 y 1
v *= fr(1) #Maximo de la funcion en el intervalo
if v < fr(u):
r.append(u)
theta = []
while len(theta) < 10000:
u = st.uniform(scale=np.pi).rvs() #Uniforme entre 0 y pi
v = st.uniform.rvs() #Uniforme entre 0 y 1
v *= 1/2 #Maximo de la funcion en el intervalo
if v < ft(u):
theta.append(u)
phi = st.uniform(scale=2*np.pi).rvs(size=10000) #Uniforme entre 0 y 2pi
hist(r, bins=20, normed=True, **formato)
x = np.linspace(0, 1, 100)
plt.plot(x, fr(x), 'g', lw=3, alpha = 0.75)
hist(theta, bins=20, normed=True, **formato)
x = np.linspace(0, np.pi, 100)
plt.plot(x, np.sin(x)/2, 'g', lw=3, alpha = 0.75)
hist(phi, bins=10, normed=True, **formato)
x = np.linspace(0, 2*np.pi, 2)
plt.plot(x, 1/(2*np.pi)*np.ones(x.shape), 'g', lw=3, alpha = 0.75)
