Guia 8, Ejercicio 13¶
import numpy as np
import scipy.optimize as opt
import scipy.stats as st
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
%matplotlib inline
np.set_printoptions(precision=3, suppress=True)
plt.rcParams['figure.figsize'] = (10, 5)
plt.rcParams['lines.linewidth'] = 3
plt.rcParams['ytick.labelsize'] = 12
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['axes.labelsize'] = 20
plt.rcParams['axes.titlesize'] = 20
plt.rcParams['legend.fontsize'] = 20
plt.rcParams['legend.numpoints'] = 1
# Cargo los datos
xdata, ydata = np.loadtxt('doble_exp.dat', skiprows=2, unpack=True)
a4, a5 = 209.69, 34.244
# Grafico
plt.subplot(111, yscale='log')
plt.errorbar(xdata, ydata, yerr=np.sqrt(ydata), fmt='o')
plt.grid(which='both')
plt.xlabel('Tiempo [s]')
plt.ylabel('Electrones detectados')
La teorÃa nos dice que el decaimiento de cada isótopo sigue una forma exponencial y que podemos modelar este fenómeno como:
$$ y(t) = a_1 + a_2 \, \exp (-t / a_4) + a_3 \, \exp (-t / a_5) $$donde $a_1$ es la radiación de fondo, $a_2$ y $a_3$ las amplitudes de los isótopos con vidas medias $a_4$ y $a_5$.
Si se considera a las vidas medias conocidas ($a_4 = 209.69 \, s$, $a_5 = 34.244 \, s$), este es un problema de ajuste lineal:
$$ \mathbf{A} \, \theta = y $$donde $\theta = (a_1, a_2, a_3)^T$, $\mathbf{A}_i = (1, \exp(-t_i / a_4), \exp(-t_i / a_5))$, e $y = (y_1 ... y_n)^T$.
Cuando $n$ (el número de mediciones) es mayor a la cantidad de parámetros, el problema se encuentra sobredeterminado, y el mejor $\theta$ en el sentido de cuadrados mÃnimos se obtiene minimizando:
$$ \chi^2 = (y - \mathbf{A} \, \theta)^T \, \mathbf{V}^{-1} \, (y - \mathbf{A} \, \theta) $$El $\theta$ mÃnimo se puede hallar analiticamente y es el siguiente:
$$ \hat{\theta} = (\mathbf{A}^T \mathbf{V}^{-1} \mathbf{A})^{-1} \mathbf{A}^T y $$$$ \mathbf{[V_\theta]}_{ij} = \text{Cov}\,(\hat{\theta}_i, \hat{\theta}_j) = (\mathbf{A}^T \mathbf{V}^{-1} \mathbf{A})^{-1}_{ij} $$donde $\mathbf{V}_{ij} = \text{Cov}(y_i, y_j)$.
En este caso, $ \text{Cov}(y_i, y_j) = \text{Var}(y_i) \, \delta_{ij} = y_i \, \delta_{ij} $, ya que el error es poissoniano.
En el caso donde no se conocen los tiempos de vida, el problema pasa a ser no lineal, por lo que el $\chi^2$ deja de ser cuadratico en los parametros (en $a_4$ y $a_5$), y la solución tiene que hallarse minimizando numericamente:
$$ \chi^2 = \sum_{i=1}^n \left(\frac{y_i - y(t_i)}{\sigma_i}\right)^2 $$def lls(A, Y, V):
'''Devuelve los parametros y su matriz de covarianza de cuadrados minimos lineales.'''
A = np.asmatrix(A)
Y = np.asmatrix(Y)
Vi = np.asmatrix(V).I
cov = (A.T * Vi * A).I
param = cov * A.T * Vi * Y
return param, cov
def modelo_lineal(t, a4, a5):
'''Devuelve la matriz del modelo.'''
y = np.empty((3, t.size))
y[0] = 1
y[1] = np.exp(-t/a4)
y[2] = np.exp(-t/a5)
return np.transpose(y)
def modelo_nolineal(t, *a):
'''Evalua la funcion y(t) = a1 + a2 exp(−t/a4) + a3 exp(−t/a5)'''
a = np.asarray(a[0]) if len(a) == 1 else np.asarray(a)
return a[0] + a[1]*np.exp(-t/a[3]) + a[2]*np.exp(-t/a[4])
Ajuste lineal¶
Del ajuste lineal se obtienen los siguientes parametros y matriz de covarianza
# Ajuste lineal
A = modelo_lineal(xdata, a4, a5)
Y = ydata[:, np.newaxis]
V = np.diag(ydata)
param, cov = lls(A, Y, V)
print('Parametros:\n', param)
print('\nMatriz de covarianza:\n', cov)
Ajuste no lineal¶
Para el ajuste no lineal, tomo como parametros iniciales a los obtenidos en el ajuste lineal y los tiempos de vida conocidos.
Los parametros encontrados son practicamente iguales a los iniciales. En cambio, son distintos los valores en la matriz de covarianza (en el bloque correspondiente los 3 parametros lineales).
# Parametros iniciales
p0 = np.ravel(param).tolist()
p0.append(a4)
p0.append(a5)
# Ajuste no lineal
popt, pcov = opt.curve_fit(modelo_nolineal, xdata, ydata, p0=p0, sigma=np.sqrt(ydata), absolute_sigma=True)
print('Ajuste no lineal 1')
print('Parametros:\n', popt)
print('\nMatriz de covarianza:\n', pcov)
# Grafico
plt.subplot(111, yscale='log')
plt.plot(xdata, A*param, label='Ajuste')
plt.errorbar(xdata, ydata, yerr=np.sqrt(ydata), fmt='o', label='Mediciones')
plt.grid(which='both')
plt.legend()
plt.xlabel('Tiempo [s]')
plt.ylabel('Electrones detectados')
Probabilidad relativa de formación¶
Para calcular la probabilidad relativa de formación, tanto el ajuste lineal como no-lineal llegan al mismo valor, pero con distinta incerteza, ya que las matrices de covarianza eran distintas (en este caso nos interesa el bloque perteneciente a $a_2$-$a_3$).
# Probabilidad relativa de formacion
def prob_rel(p0, s0):
'''
p0: parametros
s0: matriz de covarianza
'''
r = p0[1]/p0[2]
rs = r * np.sqrt(s0[1,1]/p0[1]**2 + s0[2,2]/p0[2]**2 - 2*s0[1,2]/p0[1]/p0[2])
return np.array([r, rs])
prf1 = prob_rel(np.asarray(param)[:,0], np.asarray(cov))
prf2 = prob_rel(popt, pcov)
print('Probabilidad relativa de formación\n')
print('{} | {}'.format('Tipo de ajuste'.rjust(15), 'Probabilidad relativa'))
print('-'*41)
print('{} | {:.3f} ± {:.3f}'.format('Lineal'.rjust(15), prf1[0], prf1[1]))
print('{} | {:.3f} ± {:.3f}'.format('No lineal'.rjust(15), prf2[0], prf2[1]))
Regiones de confianza¶
Cuando el problema es lineal, y los errores son gaussianos, la función $\chi^2$ tiene una distribución $\chi^2_n$ de $n$ grados de libertad, donde $n$ es la cantidad de mediciones. Esta se puede descomponer en una suma de 2 variables $\chi^2$:
$$ \begin{align} \chi^2_{n} &= \chi^2_{n-k} + \chi^2_{k} \\ (y - \mathbf{A} \, \theta)^T \, \mathbf{V}^{-1} \, (y - \mathbf{A} \, \theta) &= (y - \mathbf{A} \, \hat{\theta})^T \, \mathbf{V}^{-1} \, (y - \mathbf{A} \, \hat{\theta}) + (\hat{\theta} - \theta)^T \, \mathbf{V_\theta}^{-1} \, (\hat{\theta} - \theta) \end{align} $$donde $k$ es la cantidad de parámetros.
Entonces, las regiones de confianza para los parametros $\theta$ son elipsoides centradas en $\hat{\theta}$, definidas por la matriz $\mathbf{V_\theta}$.
Cuando los errores no son gaussianos, o el problema es no lineal, la función $\chi^2$ no tiene distribución $\chi^2_n$, pero sin embargo se considera como aproximación. Por ejemplo, en el caso de los errores poissonianos, mientras más grande el parametro de la poissoniana, más parecida es a una normal y, por lo tanto, mejor la aproximación.
Una mejor aproximación a la región de confianza se consigue dibujando el hipervolumen $\chi^2 - \chi^2_{min} < \delta$, cuyo nivel de confianza viene dado por $\int_0^\delta \chi^2_k (x) \, dx$.
Por ejemplo, para $\delta = 1$ tenemos:
x = st.chi2.cdf(1, np.arange(1, 6))
print('k | CL')
print('-'*9)
for k, p in enumerate(x, 1):
print('{} | {:.3f}'.format(k, p))
O, si queremos un CL de $68\%$ o $95\%$, los $\delta$ que queremos son:
x1 = st.chi2.ppf(.68, np.arange(1, 6))
x2 = st.chi2.ppf(.95, np.arange(1, 6))
print('k | 68% | 95%')
print('-'*17)
for k in range(x1.size):
print('{} | {:.3f} | {:.3f}'.format(k+1, x1[k], x2[k]))
En la guÃa, se muestran 2 curvas de $\delta = 1$ en el plano $a_4$-$a_5$. La más grande corresponde al contorno de la proyección del hipervolumen de 5 dimensiones sobre dicho plano. Por lo tanto, la probabilidad de que $a_4$-$a_5$ esten en esa región viene dada por la $\chi^2_5$ y es $3.7\%$. Por otro lado, la más chica corresponde a $\chi^2 - \chi^2_{min} < 1$ dejando fijos los parametros $a_1, a_2$ y $a_3$. Es decir, se modificó la distribución de probabilidad en el subespacio $a_1, a_2, a_3$ por una función delta de Dirac centrada en $\hat{a_1}, \hat{a_2}, \hat{a_3}$. Como solo hay 2 parámetros "libres", la probabilidad de que se encuentren en esa región $a_4$-$a_5$ viene dada por la $\chi^2_2$ y es $39.3\%$.
En los siguientes graficos se reproduce esas curvas, y las correspondientes a los otros parametros.
En el primero, se calcula el $\chi^2$ moviendo solo 2 parametros y dejando los otros 3 fijos. Luego se buscan las curvas de nivel $\delta = \{1, 2.279, 5.991\}$, que corresponden a CL de $0.393, 0.68$ y $0.95$, respectivamente.
Además, en la diagonal del grafico, se encuentra el $\chi^2$ moviendo solo 1 parametro y dejando el resto fijos, y se grafican lineas horizontales para $\delta = \{1, 0.989, 3.841\}$ que corresponden a CL de $0.683, 0.68$ y $0.95$, respectivamente.
Si quisiera dar intervalos de confianza de $68\%$ y $95\%$ para cada parametro por separado, serÃa en la intersección la curva azul en la diagonal con las lineas horizontales amarillas y rojas, respectivamente. Mientras que si quisiera dar intervalos de confianza conjuntos (de a 2), serÃa el area encerrada por las curvas amarillas y rojas para CL $68\%$ y $95\%$, respectivamente.
def chi2(param, xdata, ydata):
return np.sum((ydata - modelo_nolineal(xdata, param))**2/ydata)
sig = np.sqrt(np.diag(pcov))
x = np.linspace(-1, 1, 51)
a_fix = [x * sig[i] + popt[i] for i in range(popt.size)]
chi2_min = chi2(popt, xdata, ydata)
chi2_fix = dict()
for i in range(4):
for j in range(i+1, 5):
tmp = np.empty((a_fix[i].size, a_fix[j].size))
param = np.copy(popt)
for ii, ai in enumerate(a_fix[i]):
for jj, aj in enumerate(a_fix[j]):
param[[i, j]] = [ai, aj]
tmp[ii,jj] = chi2(param, xdata, ydata)
chi2_fix[(i,j)] = tmp - chi2_min
for i in range(5):
tmp = np.empty(a_fix[i].size)
param = np.copy(popt)
for ii, ai in enumerate(a_fix[i]):
param[i] = ai
tmp[ii] = chi2(param, xdata, ydata)
chi2_fix[i] = tmp - chi2_min
fig = plt.figure(figsize=(20, 10))
ax = dict()
colors = ['g', 'y', 'r']
levels1 = [1, st.chi2.ppf(0.68, 1), st.chi2.ppf(0.95, 1)]
levels2 = [1, st.chi2.ppf(0.68, 2), st.chi2.ppf(0.95, 2)]
lines = ['']*4
max_ticks = 5
for i in range(5):
for j in range(i, 5):
ax[(i,j)] = plt.subplot2grid((5, 5), (i,j))
xloc = plt.MaxNLocator(max_ticks)
yloc = plt.MaxNLocator(max_ticks)
ax[(i,j)].xaxis.set_major_locator(xloc)
ax[(i,j)].yaxis.set_major_locator(yloc)
plt.grid()
if i == j:
lines[0] = plt.plot(a_fix[i], chi2_fix[i])[0]
for c, level in enumerate(levels1):
lines[c+1] = plt.plot(a_fix[i], level * np.ones_like(a_fix[i]), color=colors[c])[0]
else:
plt.contour(a_fix[j], a_fix[i], chi2_fix[(i,j)], levels = levels2, colors = colors)
for k in range(5):
ax[(k,k)].set_ylabel('$a_' + str(k+1) + '$', rotation=0, ha='right', va='center')
ax[(0,k)].set_title('$a_' + str(k+1) + '$')
ax[(k,4)].yaxis.set_label_position('right')
ax[(k,4)].set_ylabel('$a_' + str(k+1) + '$', rotation=0, ha='left', va='center')
ax[(k,k)].set_xlim((a_fix[k][0], a_fix[k][-1]))
fig.legend(lines, ['$\chi^2_{(1)}$', '39.3%', '68%', '95%'], loc = (0.1,0.1))
En el segundo grafico (debajo), se calcula el $\chi^2$ para cada punto del plano $a_i$-$a_j$ minimizando respecto de los otros 3 parametros. Luego se buscan las curvas de nivel $\delta = \{1, 5.861, 11.070\}$, que corresponden a CL de $0.037, 0.68$ y $0.95$, respectivamente.
Además, en la diagonal del grafico, se encuentra el $\chi^2$ para 1 solo parametro minimizando respecto del resto, y se grafican lineas horizontales para $\delta = \{1, 5.861, 11.070\}$ que corresponden a CL de $0.037, 0.68$ y $0.95$, respectivamente.
Los intervalos de confianza son conjuntos a 5 parametros, y son un hipervolumen de 5 dimensiones.
def chi(paramklm, klm, param, xdata, ydata):
param[klm] = paramklm
return (ydata - modelo_nolineal(xdata, param))/np.sqrt(ydata)
sig = np.sqrt(np.diag(pcov))
x = np.linspace(-3, 3, 31)
a_free = [x * sig[i] + popt[i] for i in range(popt.size)]
chi2_min = chi2(popt, xdata, ydata)
chi2_free = dict()
for i in range(4):
for j in range(i+1, 5):
tmp = np.empty((a_free[i].size, a_free[j].size))
param = np.copy(popt)
klm = np.setdiff1d(np.arange(5), [i, j])
paramklm = param[klm]
for ii, ai in enumerate(a_free[i]):
for jj, aj in enumerate(a_free[j]):
param[[i, j]] = [ai, aj]
poptklm = opt.leastsq(chi, paramklm, args=(klm, param, xdata, ydata))[0]
newparam = np.copy(param)
newparam[klm] = poptklm
tmp[ii,jj] = chi2(newparam, xdata, ydata)
chi2_free[(i,j)] = tmp - chi2_min
for i in range(5):
tmp = np.empty(a_free[i].size)
param = np.copy(popt)
klm = np.setdiff1d(np.arange(5), i)
paramklm = param[klm]
for ii, ai in enumerate(a_free[i]):
param[i] = ai
poptklm = opt.leastsq(chi, paramklm, args=(klm, param, xdata, ydata))[0]
newparam = np.copy(param)
newparam[klm] = poptklm
tmp[ii] = chi2(newparam, xdata, ydata)
chi2_free[i] = tmp - chi2_min
fig = plt.figure(figsize=(20, 10))
ax = dict()
colors = ['g', 'y', 'r']
levels = [1, st.chi2.ppf(0.68, 5), st.chi2.ppf(0.95, 5)]
lines = ['']*4
max_ticks = 5
for i in range(5):
for j in range(i, 5):
ax[(i,j)] = plt.subplot2grid((5, 5), (i,j))
xloc = plt.MaxNLocator(max_ticks)
yloc = plt.MaxNLocator(max_ticks)
ax[(i,j)].xaxis.set_major_locator(xloc)
ax[(i,j)].yaxis.set_major_locator(yloc)
plt.grid()
if i == j:
lines[0] = plt.plot(a_free[i], chi2_free[i])[0]
for c, level in enumerate(levels):
lines[c+1] = plt.plot(a_free[i], level * np.ones_like(a_free[i]), color=colors[c])[0]
else:
plt.contour(a_free[j], a_free[i], chi2_free[(i,j)], levels = levels, colors = colors)
for k in range(5):
ax[(k,k)].set_ylabel('$a_' + str(k+1) + '$', rotation=0, ha='right', va='center')
ax[(0,k)].set_title('$a_' + str(k+1) + '$')
ax[(k,4)].yaxis.set_label_position('right')
ax[(k,4)].set_ylabel('$a_' + str(k+1) + '$', rotation=0, ha='left', va='center')
ax[(k,k)].set_xlim((a_free[k][0], a_free[k][-1]))
fig.legend(lines, ['$\chi^2_{(5)}$', '3.7%', '68%', '95%'], loc = (0.1,0.1))
Finalmente, en el último grafico, se superponen las curvas de $\delta = 1$ para cuando $\chi^2$ está fijo en el resto de los parametros (curva azul) o cuando se le permite variar a los parametros para minimizar el $\chi^2$ en cada punto del plano (curva roja).
fig = plt.figure(figsize=(20, 10))
ax = dict()
levels = [1]
lines = ['']*3
max_ticks = 5
af, bf = 7, -7
for i in range(5):
for j in range(i, 5):
ax[(i,j)] = plt.subplot2grid((5, 5), (i,j))
xloc = plt.MaxNLocator(max_ticks)
yloc = plt.MaxNLocator(max_ticks)
ax[(i,j)].xaxis.set_major_locator(xloc)
ax[(i,j)].yaxis.set_major_locator(yloc)
plt.grid()
if i == j:
lines[0] = plt.plot(a_fix[i], chi2_fix[i], color = 'b')[0]
lines[1] = plt.plot(a_free[i][af:bf], chi2_free[i][af:bf], color = 'r')[0]
for level in levels:
lines[2] = plt.plot(a_free[i][af:bf], level * np.ones_like(a_free[i][af:bf]), color='black')[0]
else:
plt.contour(a_fix[j], a_fix[i], chi2_fix[(i,j)], levels = levels, colors = 'b')
plt.contour(a_free[j][af:bf], a_free[i][af:bf], chi2_free[(i,j)][af:bf, af:bf], levels = levels, colors = 'r')
for k in range(5):
ax[(k,k)].set_ylabel('$a_' + str(k+1) + '$', rotation=0, ha='right', va='center')
ax[(0,k)].set_title('$a_' + str(k+1) + '$')
ax[(k,4)].yaxis.set_label_position('right')
ax[(k,4)].set_ylabel('$a_' + str(k+1) + '$', rotation=0, ha='left', va='center')
ax[(k,k)].set_xlim((a_free[k][af], a_free[k][bf-1]))
ax[(k,k)].set_ylim((np.min(chi2_free[k][af:bf]), np.max(chi2_free[k][af:bf])))
fig.legend(lines, ['$\chi^2_{(1)} \, ó \, \chi^2_{(2)}$', '$\chi^2_{(5)}$', '$\delta = 1$'], loc = (0.1,0.1))
