NumPy es LA biblioteca para cálculo vectorial. Además de contener un nuevo tipo de dato que nos va a ser muy útil para representar vectores y matrices, nos provee de un arsenal de funciones de todo tipo.
Vamos a empezar por importar la bliblioteca numpy. La sintaxis típica de eso era import biblio as nombre:
In [1]:
import numpy as np # con eso voy a poder acceder a las funciones de numpy a través de np.función()
# en ejemplo
print('El numero e = {:.3}'.format(np.e))
print('O el numero Pi = ', np.pi)
# y podemos calcular senos y cosenos entre otras cosas
print(np.sin(np.pi)) # casi cero! guarda con los floats!
Todo eso está muy bien, pero lo importante de numpy son los arrays numéricos. Los arrays numéricos nos van a servir para representar vectores (el objeto matemático) o columnas/tablas de datos (el objeto oriyinezco o de laboratorio).
La idea es que es parecido a una lista: son muchos números juntos en la misma variable y están indexados (los puedo llamar de a uno dando la posición dentro de la variable). La gran diferencia con las listas de Python es que los arrays de numpy operan de la forma que todos queremos:
Veamos ejemplos usando la función array para crear arrays básicos.
In [2]:
a = np.array([1, 2, 3, 4]) # array toma como argumento un vector_like (lista, tupla, otro array...vector_like)
b = np.array([5, 6, 7, 8])
print(type(a)) #Tipos de dato de a
print(a + b) # vector suma
print(a * b) # acá multiplicó
# en el caso de las listas esto era muy molesto
l1 = [1, 2, 3, 4]
l2 = [5, 6, 7, 8]
print(l1 + l2) # sumar concatena
# print(l1 * l2) # esto ni siquiera se puede hacer!
Y al igual que con las listas, uno puede acceder a elementos específicos de un array:
In [3]:
print(a[0], a[1], a[2], a[3]) # son 4 elementos, los indices van del 0 al 3
# y más o menos vale todo lo que valía con listas
print(b[-1]) # agarro al último elemento de b
print(b[1:]) # desde el primero hasta el 3 (no incluido el final, nunca se incluye)
Para facilitar la vida del usuario numpy viene con un montón de rutinas de creación de arrays típicos. En particular, matrices típicas como las identidades o de todos elementos iguales a 1 o 0 y arrays con cierta cantidad de elementos entre dos números (muy útil para crear dominios para gráficos).
Veamos ejemplos de esos:
In [4]:
# equiespaciados
equilin = np.linspace(0, 1, 10) # 10 número equiespaciados linealmente entre 0 y 1
print('Equiespaciado lineal:', equilin)
arange = np.arange(0, 1, 1./10) # como el range de Python pero puede venir con un paso en coma flotante
print('Como el range de las listas:', arange)
identidad = np.identity(3)
print('Identidad de 3x3:', identidad)
print()
#otros para que prueben ustedes
ceros = np.zeros((4, 4)) # todos ceros, matriz de 4x4
unos = np.ones((2,3)) # todos unos, matriz de 2x3
ojos = np.eye(5, k=0) # unos en la diagonal, como identidad
ojos2 = np.eye(5, k=2) # qué pasó acá?
print(ceros)
print()
print(unos)
print()
print(ojos)
print()
print(ojos2)
Y antes de seguir, algo que siempre puede ser útil: los arrays tienen ciertas propiedades como su shape (de cuánto por cuánto) y el dtype (qué tipo de cosas tiene adentro). Podemos acceder a estos datos de la siguiente manera:
In [5]:
x = np.linspace(0, 10, 1000) # ese array tiene 1000 elementos, andar printeando es poco práctico!
print(x.dtype) # array.dtype nos dice que tipo de elementos tiene el array
ceros = np.zeros((100, 100)) # matriz de 100x100
print(ceros.shape) # array.shape nos dice cuántas filas y columnas tiene el array
# prueben que pasa cuando le piden el shape a un array con una sola fila o columna como el x
Para construir un vector (o array) más complejo podemos usar las funciones r_ y c_ (busquen la documentación de c_, pero r_ mostramos un ejemplo)
In [6]:
print("Vector concatenado: ", np.r_[np.arange(-1, 1, 0.3), 2, 3, 4, np.linspace(5, 6, 10)])
#Para concatenar tres vectores podemos usarlo
a = np.linspace(0, 5, 5)
b = np.linspace(6, 11, 5)
c = np.linspace(12, 17, 5)
print("Matriz concatenada: ", np.r_['0,2', a, b, c].T)
In [7]:
# Realicen el ejercicio 1
In [8]:
# Realicen el ejercicio 2
numpy es además la librería encargada de las operaciones de álgebra lineal. Ya de por si la suma y producto por escalar de arrays de dos dimensiones son operaciones matriciales. Leanse la documentación de numpy.linalg (o revisen nuestro apartado de álgebra lineal del taller numérico); la idea es que está hecho para ser MATLAB compatible.
Hacer gráficos es lo primero que aprendimos a hacer en Origin, así que es lo primero que vamos a aprender para reemplazarlo. Van a ver que no es nada complicado.
Primero, debemos importar las bibliotecas necesarias para graficar, numpy por si no la teníamos, y de la biblioteca matplotlib (que tiene infinitas funciones y posibilidades), solamente pyplot de donde sacaremos las funciones que necesitaremos para graficar.
In [9]:
from matplotlib import pyplot as plt
# muestra los gráficos en el mismo notebook
%matplotlib inline
Importemos los datos que tenemos en el archivo datos.csv. Este archivo tiene tres columnas, x, y, errores en y. Separemoslos en tres variables
In [10]:
# Con esto cargamos los datos
# Si vemos los archivos son dos columnas
# que vamos a guardar como x, y
data = np.loadtxt("doble_exp.dat")
x = data[:,0]
y = data[:,1]
error_y = data[:,2]
Para tener un poco de intuición de los datos, sin el error en y, primeros grafiquemoslos con plt.plot
In [11]:
# Ploteamos. Los datos como puntos de color rojo
plt.plot(x, y, 'ro', label = 'Datos')
# Detalles del gráfico
plt.grid(True) # Para que quede en hoja cuadriculada
plt.title('Grafico ejemplo')
plt.xlabel('Valores en x')
plt.ylabel('Valores en y')
plt.legend(loc = 'best')
plt.show() # si no usaron %matplotlib inline, esto abre una ventanita con el gráfico
Notemos que dentro de la función plot pusimos como parámetros ro que significa que el color que queremos para la curva sea azul, y que el trazo sea una línea continua. Esto es customizable, pueden probar letras de otros colores (g, r, y, k) o bien otros trazos.
Alteremos los ejes para que sean logaritmicos (porque podemos, aunque es un proceso standard cuando tenés datos "exponenciales") y los puntos de color verde
In [12]:
# Ploteamos
plt.plot(x, y, 'go', label = 'Modelo')
# Detalles del gráfico
plt.grid(True) # Para que quede en hoja cuadriculada
plt.yscale('log')
plt.xscale('log')
plt.title('Grafico ejemplo')
plt.xlabel('Valores en x')
plt.ylabel('Valores en y')
plt.legend(loc = 'best')
plt.show()
Así aparece una recta a ajustar, con el proceso de linealizado. Esto ya lo deberían tener más o menos claro de cursar laboratorios, pero es algo muy importante; los algoritmos lineales de cuadrados mínimo (sea por la formula normal o por descendiente del gradiente) siempre convergen. No les voy a contar toda la historia porque ya lo van a ver, pero por ahora creame
Volviendo sobre la cosas para customizar los gráficos, pasense por el repositorio nuestro pero especialmente por la documentation de matplotlib
Finalmente, para graficar los errores además de los datos tenemos la función plt.errorbar (vean la documentación y entiendan por qué se usa como la usamos acá)
In [13]:
plt.errorbar(x, y, yerr=error_y, fmt='go', label="Datos")
# Detalles del gráfico
plt.grid(True) # Para que quede en hoja cuadriculada
plt.yscale('log')
plt.xscale('log')
plt.title('Grafico ejemplo')
plt.xlabel('Valores en x')
plt.ylabel('Valores en y')
plt.legend(loc = 'best')
plt.show()
Si necesitamos guardar cambios de los datos, por ejemplo con el logaritmo aplicado, podemos usar la función np.savetxt(que gaurdar en formato csv)
In [14]:
# Guardamos y cargamos, a modo de ejemplo
#vean que crea un archivo Datos_taller.txt
np.savetxt('Datos_taller.txt', np.log(data[:,:-1]), delimiter = '\t')
Data = np.loadtxt('Datos_taller.txt', delimiter = '\t')
#plt.plot(Data[:,0],Data[:,1], 'r.') # Veamos que son los mismos datos, pero linealizados!
plt.errorbar(Data[:,0],Data[:,1], yerr=1/y*error_y, fmt='ro')
plt.show()
Obtenidos los datos y con un poco de inteligencia adquirida de ellos, queremos efectuar un ajuste. Para eso importamos la biblioteca con la función que usaremos, que aplica cuadrados mínimos para obtener los coeficientes.
In [15]:
from scipy.optimize import curve_fit
El algoritmo de cuadrados mínimos necesita la función con la que queremos ajustar (que como ya linealizamos, es lineal), que vamos a definir como como función lambda, dominio, los datos, un vector con los valores iniciales de los parámetros desde donde debe comenzar a iterar. Los parámetros iniciales son importantes para el ajuste correcto (aún para el caso lineal, aunque acá el método hace un guess correcto).
La función nos devolverá 2 cosas. Primero, los parámetros optimizados por este algoritmo, ordenados como los pusimos en la función lambda cuando la definimos, que lo guardamos en el vector popt. Por otro lado nos dará la matriz de covarianza (recuerden, que tiene en su diagonal los $\sigma^2$ de cada parámetro).
In [16]:
f = lambda x, A, B: A * x + B
# Ajustamos, pero con las funciónes logaritmicas. Usamos propagación de errores
popt, pcov = curve_fit(f, np.log(x), np.log(y), sigma = 1/y * error_y,
absolute_sigma=True)
print(popt)
print(pcov)
Listo, ahora chequeamos con un gráfico que haya ajustado
In [17]:
t = np.linspace(min(x), max(x), 1000) #Las funciones "viejas" de python siguen funcionando!
t = np.log(t) #Aplico logaritmo
plt.plot(np.log(x), np.log(y), 'ro', label = 'Datos')
plt.plot(t, f(t, *popt), 'g-', label = 'Ajuste') #grafico la función
# Detalles del gráfico
plt.grid(True)
plt.title('Grafico ejemplo')
plt.xlabel('Valores en x')
plt.ylabel('Valores en y')
plt.legend(loc = 'best')
plt.show()
¿Se les ocurre otra forma de obtener el ajuste que no sea por curve_fit? Revisen la librería optimize, ahí pueden encontrar la respuesta (y con lo que van a ver en las clases)
In [18]:
# Realicen el ejercicio 3
Hasta ahora hicimos un análisis de datos adquiridos en el laboratorio, pero no hablamos de datos aleatorios (lo que manejamos cuando hablamos de estadísticas!)
La computadora tiene internamente un mecanismo de creación de números pseudoaleatorios (usando un algoritmo bastante elegante) que permite obtener muestras de una variable con distribución uniforme "continua" (¿me explican por qué las comillas?). Sobre eso, Python crea un montón de librerías para muestrear variables en otras distibuciones
Y luego, vamos a crear histogramas, es decir un gráfico donde dividimos una serie de datos en rangos y contamos cuántos de nuestros datos caen en cada rango. A esos rangos se los llama bins.
La función que vamos a usar es plt.hist toma como argumentos un array de números, en cuántos bins queremos dividir a nuestro eje x y algunas otras opciones de color como constante de normalización y color de las barras.
Hagamos un histograma simple de un set gaussiano. Para eso, creemos datos alrededor de algún valor medio usando randn de NumPy:
In [19]:
mu, sigma = 100, 15 # mu es mi valor medio, sigma la desviación
x = mu + sigma * np.random.randn(10000) # Ya deben saber, que z = (x - mu)/s es una N(0,1)
n, bins, patches = plt.hist(x, bins=50, normed=1, facecolor='green',
edgecolor='black', alpha=0.75)
# en la variable n se encuentran los datos del histograma
# bins es un vector con los bordes de los rangos de datos
# patches no nos interesa en general
Y ya que estamos, para mostrar cómo afecta la elección de bins, graficamos dos histogramas uno arriba del otro.
In [20]:
n, bins, patches = plt.hist(x, bins=100, normed=1, facecolor='green', alpha=0.75)
n, bins, patches = plt.hist(x, bins=10, normed=1, facecolor='red', alpha=0.75)
Los bins los podemos construir como un array, con linspace, pero tenemos que recordar que estamos dando el borde derecho de los bins
In [21]:
bins = np.arange(0, 200, 5)
plt.hist(x, bins = bins, normed=1, facecolor='green',
edgecolor='black', alpha=0.75);
No solo podemos muestrear variables aleatorias con distribución normal, numpy.random tiene una gran selección de distribuciones. Aún así en la librería scipy.stats tenemos aún más distribuciones y funciones estadísticas, pero requieren otra forma de uso. Veamos como muestrar una distribución $\chi^2$ con 3 grados de libertad
In [22]:
import scipy as sp
from scipy import stats #Tengo que importar la sublibrería
X = sp.stats.chi2(df=3)
t = np.linspace(0, 10, 1000)
plt.plot(t, X.pdf(t), 'r-') #Puedo usar la PDF de la distribución así
Out[22]:
La variable X que construimos con scipy.stats tiene propiedades útiles para la inferencia estadística. Es más veamos su tipo y además usemos la función help para ver que contiene
In [23]:
print(type(X))
#print(help(X)) #Hagan esto si están en un entorno no tan cool como Jupyter
El tipo rv_frozen se llama así porque es una distribución congelada en sus parámetros. Podríamos construir una distribución para un caso puntual a partir de la función scipy.stats.chi2, variando df en cada caso. Pero para mejorar la lectura, podemos considerar a X como una variable aleatoria, donde tenemos conocimiento de su distribución, de sus momentos y podemos muestrar de ella.
Veamos primero los primeros momentos y la capacidad de integrar funciones arbitrarias
In [24]:
print(X.mean(), X.moment(1))
print(X.var(), X.moment(2)) #Calculen si está bien la varianza y el momento no centrado
print(X.expect(lambda x: x**3), X.moment(3)) #Integral entre 0 e infinito de x**3.
Mientras si queremos muestrar de la distribución tenemos la función rvs
In [25]:
A = X.rvs(10000)
plt.hist(A, bins=50, normed=1, facecolor='green',
edgecolor='black', alpha=0.75);
Muestren de la distribución exponencial con scipy.stats (que pueden ver acá) y tomen 10000 muestras de esta distribución. ¿Se les ocurre una forma de encontrar una exponencial desde una uniforme? Verifiquen muestreando 1000 datos de la uniforme y aplicando la transformación que realmente es una exponencial
In [26]:
# Realicen el ejercicio 4
Con esto último ya tienen todas las herramientas para lidiar con las guías y el parcial, pero antes de terminar vamos a ver como agregarle barras de error a los histogramas.
Para hacer un gráfico de histograma con errores tenemos un par de opciones. La más intuitiva es usar las herramientas pandas, que permiten crear a partir de un histograma computado con numpy.histogram un gráfico de barras. La función, que pueden usar y adaptar para esto es
In [27]:
def plot_histogram(data, bins, ticks=5, xlabel='X', ylabel='Histograma', density=True, ecolor='', ax = None):
N = len(data)
num_bins = len(bins)
xmin = min(bins)
xmax = max(bins)
if ax is None:
ax = plt.gca()
hist, _ = np.histogram(data, bins=bins, density=density)
normed = N if density else 1
error = np.sqrt(hist / normed) # por qué es así el error?
# Se le puede pasar un yerr y xerr, como a errorbar
# ecolor es el color de los errores y capsize es el tamaño del fin de barra de error
plt.bar(x=bins[:-1], height=hist, width=np.diff(bins), yerr=error, edgecolor='black', ecolor="red", capsize=2)
plt.xticks(np.linspace(0, num_bins - 2, ticks), np.linspace(xmin, xmax, ticks), rotation='horizontal');
plt.xlabel(xlabel);
return ax
plt.figure(figsize=(7 * 1.78, 7))
ax = plot_histogram(A, np.linspace(0, 25, 50))
plt.grid(True, axis='y')
El método de Monte Carlo en estadística pertenece a una familia de métodos de remuestro, que permite obtener, de muestras pequeñas de datos, estimadores con sus valor estimado y varianza. Veamos en qué consiste este método.
Primero volvamos a cargar los datos de double_exp.dat
In [28]:
# Volvemos a cargar los datos y los guardamos en x, y
data = np.loadtxt("doble_exp.dat")
x = data[:,0]
y = data[:,1]
error_y = data[:,2]
El método consiste en simular un nuevo conjunto de mediciones usando las propiedades del error, que en este caso asumimos gaussiano (podría ser de otra propiedad, pero por el teorema central del límite casi siempre será una distribución normal).
In [29]:
#Creo al azar datos dentro del error de y con distribución gaussiana
yp = y + np.random.randn(y.shape[0]) * error_y
plt.plot(x, yp, 'ro')
plt.yscale('log')
plt.grid()
plt.errorbar(x, y, yerr=error_y, fmt='b*');
In [30]:
#Creo al azar datos dentro del error de y con distribución gaussiana
M = []
logx = np.log(x)
f = lambda x, A, B: A * x + B
for i in range(10000):
#Vuelvo a computar y
yp = y + np.random.randn(y.shape[0]) * error_y
logy = np.log(yp)
#Dado que calculo el logaritmo, debo sacarme los valores NaN (que son logaritmos negativos)
valid_idx = ~np.isnan(logy)
p, cov = sp.optimize.curve_fit(f, logx[valid_idx], logy[valid_idx])
M.append(p)
M = np.array(M)
#Valor medio y covarianza
p = np.mean(M, axis=0)
cov = np.cov(M.T)
#Imprimamos
print(p)
print(cov)
t = np.linspace(min(x), max(y), 1000)
plt.grid()
plt.errorbar(np.log(x), np.log(y), yerr=1/y * error_y, fmt='go')
plt.plot(np.log(t), f(np.log(t), *p), 'r-');
Si revisan las cuentas de curve_fit les da un resultado sumamente parecido, lo que nos da un poco más de confianza en el método (o curve_fit hace algo parecido, fijense que scipy tiene a disposición el código fuente de forma bastante a mano).
Finalmente hagamos el histograma de ambos estimadores en 2D. Para eso usamos numpy.histogram2d que nos devuelve las cuentas, las coordenadas en x, y. Lean lo que dice histogram2d para entenderlo (y por qué usamos los bordes xedges e yedges con un elemento menos)
In [31]:
counts, xedges, yedges = np.histogram2d(M[:,0], M[:,1], bins=100, normed=True);
plt.xlim((min(M[:,0]), max(M[:,0])))
plt.ylim((min(M[:,1]), max(M[:,1])))
plt.xlabel("A")
plt.ylabel("B")
plt.contour(xedges[:-1], yedges[:-1], counts);
Ahora la pregunta que deben hacerse, ¿están o no correlacionado los parámetros? ¿Tiene sentido?.
Prueben ejecutar el ajuste no linealizado con Monte Carlo. Es decir en vez de ajustar una lineal, deben ajustar $$f(x|A,B,C,D, E) = A + B \exp(C x) + D \exp(E x)$$ Elija como parámetros iniciales $A=10$, $B=130$, $C=-0.001$, $D=960$ y $E=-0.02$
In [32]:
#Ejercicio 5
Bueno, con esto terminamos las sesiones del taller de Python. Esperemos que les haya servido, y que les sirva para avanzar en la materia (y en la vida profesional). A nosotros nos sirve y trabajamos con esto todo el tiempo
Como siempre, tenemos en el reposotorio de Github (https://github.com/fifabsas/talleresfifabsas) este material y más cosas que vamos subiendo con el tiempo.
Páginas para leer
http://pybonacci.org/2012/06/07/algebra-lineal-en-python-con-numpy-i-operaciones-basicas/
http://relopezbriega.github.io/blog/2015/06/14/algebra-lineal-con-python/
http://pendientedemigracion.ucm.es/info/aocg/python/modulos_cientificos/numpy/index.html
Documentación
https://docs.python.org/2/library/math.html
http://docs.scipy.org/doc/numpy/reference/routines.linalg.html
http://matplotlib.org/api/pyplot_api.html
http://pandas.pydata.org/pandas-docs/stable/ (no la vimos, pero estos son los dataframes)
Todo esto es posible gracias al aporte de mucha gente.