import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
%matplotlib inline
a) Implementar un programa que genere el resultado de una serie de n experimentos de Bernoulli independientes y cuente el numero de exitos obtenidos.¶
Ya existen esas funciones en Scipy. Ver scipy.stats.bernoulli y scipy.stats.binom. De todas maneras, hago una implementación simple:
def Bernoulli(p):
x = stats.uniform.rvs()
if x < p:
return 1
else:
return 0
def Binomial(n, p):
x = 0
for k in range(n):
if Bernoulli(p):
x += 1
return x
b) Inciden n = 15 fotones sobre un detector con eficiencia p = 0.75. Cuantos son detectados en un experimento en particular? Repita el experimento 1000 veces, realice un histograma y comparelo con la distribucion de probabilidad teorica.¶
La distribución teórica de este experimento es una binomial con parametros $n = 15$ y $p = 0.75$.
Para la construcción del histograma, se tomaron 1000 datos aleatorios de esa distribución. El histograma se encuentra normalizado. Las barras de error del histograma se explican en el punto f).
n, p = 15, 0.75 # n: cantidad de fotones. p: eficiencia del detector
rep = 1000 # rep: repeticiones del experimento
binom = stats.binom(n, p) # Freezeo la binomial con parametros n y p.
x1 = binom.rvs(size=rep) # Tomo rep valores aleatorios de la binomial.
y1, bins = np.histogram(x1, bins=np.arange(0,n+2))
y1std = np.sqrt(rep*y1/rep*(1-y1/rep))
y1 = y1 / rep
y1std = y1std / rep
bincenters = bins[0:-1]
x2 = np.arange(0, n+1)
y2 = binom.pmf(x2) # Densidad de probabilidad de la binomial.
fig, ax = plt.subplots(1,1)
ax.plot(x2, y2, lw=5, alpha=0.6)
ax.bar(bincenters, y1, color='r', yerr=y1std, align='center')
c) Considere una fuente de intensidad media $I = 15 \frac{fot}{s}$. Simule el número de fotones en $\Delta t = 1 s$ dividiendo el intervalo en $m$ ($= 1000$) subintervalos $dt$. Aproxime la probabilidad de emitir un foton en $dt$ como $Idt$ y desprecie la probabilidad de emitir más de 1 foton en el intervalo. Justifique esta hipotesis. Simule el numero de fotones emitidos en cada $dt$ a partir del programa implementado en el punto a), y sumando sobre todos los $dt$ calcule el numero total de fotones emitidos en $\Delta t$. Repita el experimento 1000 veces, construya un histograma y comparelo con la distribución teorica.¶
La probabilidad de emisión de $k$ fotones en un tiempo $t$ de una fuente de intesidad media $I$ sigue una distribucíon poissoniana con parametro $\lambda = It$:
$$ Poisson(\lambda, k) = \frac{\lambda^k}{k!} e^{-\lambda} $$La probabilidad de emitir más de un fotón en el intervalo $dt = \Delta t / m$ se puede calcular como:
$$ P (k > 1) = \sum_{k=2}^{\infty} Poisson(\lambda = I dt, k) = 1 - Poisson(I dt, 0) - Poisson(I dt, 1) \approx 10^{-4} $$Al ser esta una probabilidad muy baja, es una buena aproximación pensar que en un intervalo $dt$ se puede emitir como mucho 1 foton, y por lo tanto hay 2 eventos: que se emita un foton o que no se emita. Es decir, es un experimento de Bernoulli.
Por otro lado, la probabilidad de que emita 1 foton es $Poisson(\lambda = Idt, 1)$. Como $\lambda = I\,dt$ es chiquito, podemos aproximar $Poisson(\lambda, 1) = \frac{\lambda^1}{1!} e^{-\lambda} = \lambda \, [1 + O(\lambda)] \approx \lambda$.
Entonces, la cantidad de fotones emitidos en $\Delta t = 1 s$ se aproxima por la suma 1000 experimentos de Bernoulli con probabilidad $p = I\,dt$.
Abajo se grafica el histograma generado y se superpone con la distribución teorica.
Int = 15 # 15 fotones / segundo
T = 1 # 1 segundo
N = 1000 # Numero de intervalos
dT = T / N # Intervalos de tiempo
rep = 1000
bernoulli = stats.bernoulli(Int*dT) # Bernoulli con parametro Int*dT
poisson = stats.poisson(Int) # Poisson con parametro Int.
x3 = np.empty(rep)
for i in range(rep):
x3[i] = np.sum(bernoulli.rvs(size=N))
y3, bins = np.histogram(x3, bins=np.arange(0,30))
y3norm = np.sum(y3)
y3 = y3 / y3norm
y3std = np.sqrt(y3norm*y3*(1-y3)) / y3norm
x4 = np.arange(0, np.max(x3))
y4 = poisson.pmf(x4)
fig, ax = plt.subplots(1,1)
ax.plot(x4, y4, lw=5, alpha=0.6)
ax.bar(bins[0:-1], y3, color='r', yerr=y3std, align='center')
fotonesMultiples = stats.poisson(Int*dT).sf(1) # La probabilidad de emitir mas de 1 foton en dT
print('Probabilidad de emitir más de 1 foton en dt = ', fotonesMultiples)
d) Dado el resultado de cada uno de los 1000 experimentos del punto anterior, calcule el numero de fotones detectados a partir del programa implementado en el punto b). Realice el histograma correspondiente y comparelo con la distribución teórica.¶
Para cada resultado ($n$ fotones emitidos) de los 1000 experimentos, se extrae un numero aleatorio de una distribucion binomial con parametros $n$ y $p = 0.75$. Se genera un histograma que corresponde a la combinacion de emitir y detectar un foton.
La distribución teorica esperada en este caso es otra poissoniana con parametro $\lambda = 0.75 \, I \Delta t$. En el punto e) se explica la composicion de una poissoniana con una binomial.
poisson = stats.poisson(Int*p)
x5 = stats.binom.rvs(x3.astype(int), p)
x6 = np.arange(0, np.max(x5))
y6 = poisson.pmf(x6)
y5, bins = np.histogram(x5, bins=np.arange(0,30))
y5norm = np.sum(y5)
y5 = y5 / y5norm
y5std = np.sqrt(y5norm*y5*(1-y5)) / y5norm
fig, ax = plt.subplots(1,1)
ax.plot(x6, y6, lw=5)
ax.bar(bins[0:-1], y5, color='r', yerr=y5std, align='center')
e) Se podría haber considerado la emisión y detección en forma conjunta, suponiendo una probabilidad efectiva para este proceso. Realice la simulación de este modo, grafique los histogramas y sus correspondientes distribuciones teóricas. Discuta el procedimiento y los resultados a la luz de la composición de un proceso de Poisson con uno de Bernoulli.¶
Siendo n la cantidad de fotones emitida, para cada n fijo, la probabilidad de detectar k fotones es una binomial con parametros n y p, siendo p la eficiencia del detector:
$$ Binom(n,p,k) = \binom{n}{k} p^k (1-p)^{n-k} = \frac{n!}{k!(n-k)!} p^k (1-p)^{n-k} $$Pero la cantidad de fotones emitidos no esta fija, sigue una distribución poissoniana con parametro $\lambda$:
$$ Poisson(\lambda, n) = \frac{\lambda^n e^{-\lambda}}{n!} $$Para ver la probabilidad efectiva del proceso de emisión y detección, combinamos ambas distribuciones sumando la binomial para cada n con un peso poissoniano:
$$ P_{efectiva} (\lambda, p, k) = \sum_{n=0}^\infty Binom(n, p, k) Poisson(\lambda, n) = \sum_{n=0}^\infty \frac{n!}{k!(n-k)!} p^k (1-p)^{n-k} \frac{\lambda^n e^{-\lambda}}{n!} = $$$$ = \frac{p^k}{k!} e^{-\lambda} \sum_{n=0}^\infty \frac{(1-p)^{n-k} \lambda^n}{(n-k)!} = \frac{p^k}{k!} e^{-\lambda} \sum_{m=0}^\infty \frac{(1-p)^{m} \lambda^{m+k}}{m!} = $$$$ = \frac{(\lambda p)^k}{k!} e^{-\lambda} \sum_{m=0}^\infty \frac{(1-p)^{m} \lambda^m}{m!} = \frac{(\lambda p)^k}{k!} e^{-\lambda} e^{(1-p) \lambda} = $$$$ = \frac{(\lambda p)^k}{k!} e^{-\lambda p} = Poisson(\lambda p, k) $$Se llega a que la probabilidad efectiva es otra poissoniana con parametro $\lambda p$.
poisson = stats.poisson(Int*p)
x7 = poisson.rvs(1000)
x8 = np.arange(0, 30)
y8 = poisson.pmf(x8)
y7, bins = np.histogram(x7, bins=np.arange(0,30))
y7norm = np.sum(y7)
y7 = y7 / y7norm
y7std = np.sqrt(y7norm*y7*(1-y7)) / y7norm
fig, ax = plt.subplots(1,1)
ax.plot(x8, y8, lw=5)
ax.bar(bins[0:-1], y7, color='r', yerr=y7std, align='center')
f) Que distribucion espera para el numero de datos en una determinada clase de histograma y por que? Use la desviacion estandar como una estimacion de la incerteza de dicha variable y grafique barras de error sobre los histogramas realizados.¶
El número de datos en cada clase 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. Si se conoce la distribución que siguen los datos con los que se construye el histograma, el cálculo de la probabilidad $p$ se puede hacer explícitamente como la probabilidad del intervalo que cubre la categoría en cuestión. Si en cambio no se conoce la distribución, una buena aproximación de la probabilidad es la $k/N$, es decir, la cantidad de elementos en esa categoría dividido la cantidad total de elementos del histograma.
Por lo tanto, la desviación estándar de la distribución binomial, $\sigma = \sqrt{N \, p \, (1-p)}$ es una buena estimación del error de cada categoría del histograma.