Taller de Python Capítulo 2: Numérico y Laboratorio¶
Nuestra motivación para hoy¶
De la última clase, tenemos herramientas básicas de Python que tienen muchisimas aplicaciones. Hoy vamos a ver como se puede usar Python para analizar datos de una caida libre. Digamos que en su laboratorio adquirieron los datos de la posición de una pelotita en distintos tiempos, luego de haberla soltado a una altura de 100 metros. Ahora, quieren cargar estos datos a un programa de Python, visualizarlos, y obtener el valor de la gravedad $g$.
Vamos a ver algunas herramientas que nos permitirán llevar esto a cabo.
Bibliotecas¶
Ya hemos visto la vez anterior que Python tiene varias funciones que vienen "de fábrica" como help(), print() así también como operaciones básicas entre números como sumar, restar, etc; también vimos que nosotros podemos crear nustras propias funciones para que hagan lo que necesitemos usando la palabra clave def nombre_funcion. Si uno necesita siempre realizar las mismas operaciones, reutilizará la misma función en todos sus códigos.
Por ejemplo: Supongamos que uno quiere calcular cos(x). Python no viene por defecto con esa operación, uno debería crear un algoritmo (es decir, una serie de acciones) que calcule el valor del cos de x (cosa que puede ser no trivial), pero es obvio que alguien ya lo hizo antes, alguien ya pensó el algoritmo, lo escribió, y lo utiliza diariamente, si esta persona subió su código a internet, todos podemos aprovechar y utilizarlo sin preocuparnos en cómo hizo esta persona!! Solamente hay que decirle a Python donde es que está guardado esta función. Esta posibilidad de usar algoritmos de otros es fundamental en la programación, porque es lo que permite que nuestro problema se limite solamente a entender cómo llamar a estos algoritmos ya pensados y no tener que pensarlos cada vez.
Vamos entonces a decirle a Python que, además de sus operaciones de fábrica, queremos ampliar nuestro abanico de operaciones matemáticas en todas las opciones que aparecen dentro de la biblioteca "math" (una biblioteca es, tal como sugiere el nombre, un archivo con un montón de funciones, objetos y demás).
import math # Importo a mi programa todo lo que este contenido dentro del archivo "math"
Con esta linea Python entiende que queremos que traiga todo lo que está dentro del archivo "math"
OJO: No todas las bibliotecas vienen instaladas por defecto, si queremos usar una muy rara es probable que tengamos que instalarla nosotros. Las que vamos a utilizar en el curso ya vienen instaladas por Google en Colaboratory o fueron instaladas automaticamente por Anaconda en su computadora.¶
Macanudo, ahora que Python trajo esa biblioteca, nosotros podemos acceder a su contenido usando la sintaxis math.lo_que_haya_dentro_de_math por ejemplo, math.cos()
math.cos(0)
1.0
Hay una infinidad de bibliotecas o librerías dando vueltas por internet. Muchas veces el problema que queremos solucionar se reduce simplemente a encontrar la librería que tenga la funcionalidad correcta.
Ahora continuaremos usando una muy utilizada en nuestro ámbito, la querídisima NumPy.
NumPy : vectores, matrices y tablas de datos¶
NumPy (NUMeric PYthon) es LA biblioteca para operar sobre vectores de números. 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:
import numpy as np # con eso voy a poder acceder a las funciones de numpy a través de np.función()
# Por ejemplo
print(f"El numero e = {np.e}")
print(f"O el numero Pi = {np.pi}")
El numero e = 2.718281828459045 O el numero Pi = 3.141592653589793
# Podemos calcular senos y cosenos de numeros igual que con math!
print(np.sin(np.pi)) # casi cero! guarda con los floats!
1.2246467991473532e-16
Todo eso está muy bien, pero lo importante de numpy son los arrays numéricos, que van a ser como una lista de números pero con muchos esteroides. Los arrays numéricos nos van a servir para representar vectores (el objeto matemático de "la tira de números", no el físico) o columnas/tablas de datos (el objeto oriyinezco o de laboratorio).
(Ojo, los arrays no son vectores, aunque pueden comportarse como tales en muchos sentidos)
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:
- Si sumamos o restamos dos arrays, se suman componente a componente.
- Si multiplicamos o dividimos dos arrays, se multiplican o dividen componente a componente.
Veamos ejemplos usando la función np.array para crear arrays básicos.
array_1 = np.array([1, 2, 3, 4])
array_2 = np.array([5, 6, 7, 8])
print("ARRAYS:")
print(f"{array_1} + {array_2} = {array_1 + array_2}")
print(f"{array_1} * {array_2} = {array_1*array_2}")
ARRAYS: [1 2 3 4] + [5 6 7 8] = [ 6 8 10 12] [1 2 3 4] * [5 6 7 8] = [ 5 12 21 32]
En el caso de las listas esto era muy molesto
lista_1 = [1, 2, 3, 4]
lista_2 = [5, 6, 7, 8]
print("LISTAS:")
print(f"{lista_2} + {lista_1} = {lista_1 + lista_2}") # sumar concatena
# print(lista_1 * lista_2) # esto ni siquiera se puede hacer!
LISTAS: [5, 6, 7, 8] + [1, 2, 3, 4] = [1, 2, 3, 4, 5, 6, 7, 8]
Y al igual que con las listas, uno puede acceder a elementos específicos de un array con su índice:
print(array_1[0], array_1[1], array_1[2], array_1[3])
# Y más o menos vale todo lo que valía con listas
print(array_2[-1]) # agarro al último elemento de b
1 2 3 4 8
También podemos iterar sobre ellos con for , de las mismas formas que habíamos visto para iterar listas.
# Aca itero sobre los ELEMENTOS de array_1
print("array_1:")
for num in array_1:
print(num)
print()
# Aca itero sobre los INDICES de array_2
print("array_2:")
for j in range(len(array_2)):
print(array_2[j])
array_1: 1 2 3 4 array_2: 5 6 7 8
Arrays de dos dimensiones¶
Para crear arrays de dos dimensiones (o más), podemos aplicar la función np.array sobre una lista de listas:
# Cada lista corresponde a una fila de la matriz
M = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(M)
[[1 2 3] [4 5 6] [7 8 9]]
Aca hacer M[0] dara la primer FILA de la matrix.
print(M[0])
[1 2 3]
Si queremos sacar la primer componente habra que hacer M[0, 0]
print(M[0, 0])
1
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:
# N = 11 puntos en el rango [0, 1] con 1/(N-1) de espaciamiento
equi_linspace = np.linspace(0, 1, 10)
print(f'Numero de puntos fijo: {equi_linspace}')
Numero de puntos fijo: [0. 0.11111111 0.22222222 0.33333333 0.44444444 0.55555556 0.66666667 0.77777778 0.88888889 1. ]
# Elijo el paso entre [0, 1), lo que me da
# {0, paso, 2*paso, ... n*paso} siempre que n*paso<1
equi_arange = np.arange(0, 1, 0.1)
print(f'Paso entre puntos fijo: {equi_arange}')
Paso entre puntos fijo: [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]
Otros ejemplos
# Matriz identidad de 3x3
identidad = np.identity(3)
print(f'Matriz identidad:\n{identidad}')
print()
# Matriz de todos 0 de 4x4
ceros = np.zeros((4, 4))
print(f'Matriz de ceros:\n{ceros}')
print()
# Matriz de todos 1 de 2x3
unos = np.ones((2,3))
print(f'Matriz de unos:\n{unos}')
print()
# Matriz con diagonal corrida de 5x5
ojos = np.eye(5, k=1)
print(f'Identidad corrida:\n{ojos}')
print()
# Matriz con la diagonal que yo le doy
diagonal = np.diag([1, 2, 3, 4])
print(f'Matriz con diagonal cualquiera:\n{diagonal}')
Matriz identidad: [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]] Matriz de ceros: [[0. 0. 0. 0.] [0. 0. 0. 0.] [0. 0. 0. 0.] [0. 0. 0. 0.]] Matriz de unos: [[1. 1. 1.] [1. 1. 1.]] Identidad corrida: [[0. 1. 0. 0. 0.] [0. 0. 1. 0. 0.] [0. 0. 0. 1. 0.] [0. 0. 0. 0. 1.] [0. 0. 0. 0. 0.]] Matriz con diagonal cualquiera: [[1 0 0 0] [0 2 0 0] [0 0 3 0] [0 0 0 4]]
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:
# Array de 1000 elementos
x = np.linspace(0, 10, 1000)
# ¿De qué tipos son los elementos de x?
print(x.dtype) # Ojo, sin paréntesis
float64
# Matriz llena de 0 de 40x100
ceros = np.zeros((40, 100))
# De qué tamaño es la matriz
print(ceros.shape)
(40, 100)
Ejercicio 0¶
Vamos a crear los valores de tiempo y posición que uno espera para una caida libre. Para esto, vamos a necesitar un array de valores 10 equiespaciados para el tiempo $t$, desde 0 hasta 3 segundos, y luego obtener la posición en metros mediante la ecuación:
$$y(t) = y_0 - \frac{g}{2}t^2$$Donde asumimos velocidad inicial nula, $g = 9.81 \mathrm{\frac{m}{s^2}}$ e $y_0 = 100\; \mathrm{m}$. Llamar las variables de tiempo $t$ e $y$ como t_teorico e y_teorico respectivamente.
Mas adelante compararemos estos valores con otros medidos en el laboratorio.
Recordar que podemos aplicar operaciones matemáticas a los arrays
Llevando los datos a Python¶
Muchas veces estamos acostumbrados a guardar los datos en hojas de Excel o en las Hoja de Cálculo de Drive. Python puede leer estas hojas utilizando las librerías adecuadas, pero es más sencillo, y es mejor práctica, exportar estos archivos en CSV para luego leerlos desde ahí.
Por qué CSV¶
- El formato CSV (Comma Separated Values) es abierto y un estándar a la hora de guardar datos de tipo fila o columna. Guardar los datos en Excel o Drive no está mal, pero se dificulta a la hora de pasarle los archivos a un/a colega, los CSV los puede abrir cualquiera! Inclusive, se pueden abrir Excel y Drive si se busca visualizarlo mejor.
- Con pocos datos no hay problemas de lag o crasheos, pero si son muchos datos: "Microsoft Excel dejó de responder" .$^1$
Acá una imagen de cómo descargar los archivos de Drive en CSV y de cómo se ven una vez en la compu. Es solo un archivo de texto con valores separados con comas, nada místico.
$\scriptsize\text{1. Para archivos muchos muuuchos datos inclusive los CSV se pueden quedar corto, y para estos casos existen formatos específicos que dependen de la naturaleza de los datos.}$
Una vez tenemos el archivo que descargamos podemos leerlo desde Python utilizando numpy. En este caso, la linea es
data = np.loadtxt('nombre_del_archivo.csv', delimiter=',', skiprows=1, unpack=True)
Desglosemos esto por argumentos:
El primer argumento es ubicación del archivo con la ruta respecto del archivo de Python. En este caso tenemos ambos archivos en la misma carpeta por lo que es solo el nombre.
El segundo argumento, delimiter, recibe un string que le indica a Python cómo están separados los valores en nuestro archivo. En este caso es por comas porque son archivos CSV.
El tercer argumento es más opcional, y nos permite evitar que Python lea algunas lineas del archivo. En nuestro caso pusimos 1 ya que queremos evitar que lea la primera linea, que dice "VAR 1, VAR 2".
Por qué
unpack=Truenp.loadtxt()por defecto te da los valores exactamente como están en el CSV, por lo que en este caso nos daría 2 columnas de datos con 20 filas. Utilizandounpack=Truetrasponemos los datos para pasar las columnas a filas, de modo que la primera fila (data[0]) sean los datos de VAR 1 y la segunda fila (data[1]) sean los valores de VAR 2. Acá el código utilizandounpack=Truey sin utilizarloTambién podríamos importarlo sin usar el
unpack=Truey luego usardata = np.transpose(data)
data = np.loadtxt('data_exp.csv', delimiter=',', skiprows=1)
print(f"Data sin unpack:\n {data}")
print(f'Primera columna:\n{data[0]}')
Data sin unpack: [[ 0. 1.93 1. ] [ 1.58 1.67 1. ] [ 3.16 1.3 1. ] [ 4.74 0.63 1. ] [ 6.32 1.26 1. ] [ 7.89 2.38 1. ] [ 9.47 2.91 1. ] [11.05 2.36 1. ] [12.63 3.93 1. ] [14.21 4.8 1. ] [15.79 4.83 1. ] [17.37 5.17 1. ] [18.95 6.21 1. ] [20.53 7.94 1. ] [22.11 9.12 1. ] [23.68 11.09 1. ] [25.26 12.57 1. ] [26.84 14.28 1. ] [28.42 17.1 1. ] [30. 20.56 1. ]] Primera columna: [0. 1.93 1. ]
data = np.loadtxt('data_exp.csv', delimiter=',', skiprows=1, unpack=True)
print(f"\nData con unpack:\n {data}")
print(f'Primera columna:\n{data[0]}')
Data con unpack: [[ 0. 1.58 3.16 4.74 6.32 7.89 9.47 11.05 12.63 14.21 15.79 17.37 18.95 20.53 22.11 23.68 25.26 26.84 28.42 30. ] [ 1.93 1.67 1.3 0.63 1.26 2.38 2.91 2.36 3.93 4.8 4.83 5.17 6.21 7.94 9.12 11.09 12.57 14.28 17.1 20.56] [ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. ]] Primera columna: [ 0. 1.58 3.16 4.74 6.32 7.89 9.47 11.05 12.63 14.21 15.79 17.37 18.95 20.53 22.11 23.68 25.26 26.84 28.42 30. ]
Ejercicio 1¶
Ahora vamos a cargar los datos guardados en el laboratorio, utilizando un .csv y las funcionalidades de Numpy. Guardar el tiempo y la posición en dos variables t_labo e y_labo.
Para cargar el archivo se pueden seguir los siguientes pasos:
- Descargar el archivo en este link. Pueden hacerlo con click derecho -> "Guardar como..."
- Importen el archivo data_caida_libre.csv a Python (miren el archivo para ver cuántas filas tienen que saltarse)
Graficar funciones¶
A lo largo de nuestras carreras, y de la vida, nos encontramos con información que queremos graficar, como por ejemplo la posición en el tiempo de la pelota que cae. Vamos a ver que con la ayuda de la biblioteca matplotlib esto no es nada complicado.
Esta librería es muy grande, por lo que tiene divididas sus funciones en distintos módulos. Para lo que vamos a ver en esta segunda parte nos alcanza con la sublibrería pyplot. Una de las formas de importar esta librería es:
import matplotlib.pyplot as plt
Podemos usar esta librería en conjunto con numpy para graficar las funciones que querramos.
Veamos un ejemplo, grafiquemos la función $f(x) = \sin(x)\sin(20x)$ entre $0$ y $2\pi$ (no tienen por qué entender la función, es solo el logo de Arctic Monkeys)
import numpy as np
import matplotlib.pyplot as plt
# Definimos la función
def f(x):
return np.sin(x)*np.sin(20*x)
# Definimos el dominio y la imagen
x = np.linspace(0, 2*np.pi, 500)
y = f(x)
# Graficamos la función
plt.plot(x, y)
plt.show()
Desglosemos el código que acabamos de escribir:
- Primero escribimos la función que queremos graficar.
- Luego, definimos nuestro dominio
xy nuestra imageny. - Y por último usamos la función de
plotdematplotlib.pyplotpara graficaryen función dex.
Aclaración 1: plt.show(): Si bien esta linea no es siempre necesaria en los notebooks (como google Colab), su propósito es "cerrar" la figura y mostrarla.
Aclaración 2: y no es una función. Sino que es un array que contiene los números que resultan de aplicarle f a cada uno de los valores de x. Lo que vemos arriba parece ser una función continua pero en realidad es un conjunto de punto unidos con una línea. Para visualizar mejor esto podemos hacer un gráfico que muestre explícitamente los puntos que forman el gráfico:
plt.plot(x, y, '.-')
plt.show()
Qué hace plt.plot()?¶
Esta función que parece mágica lo único que hace es tomar los vectores que le damos como inputs y graficar un punto en cada par de coordenadas $(x,y)$. Además, por defecto une estos puntos con lineas rectas. Por lo tanto, para que la función pueda hacer su trabajo, es necesario que los array x e y tengan la misma longitud. Para dejar claro este concepto veamos un ejemplo con pocos puntos escritos a mano:
coords_en_x = np.array([0, 2, 6, 8, 10])
coords_en_y = np.array([2, 4, 6, 14, 14])
plt.plot(coords_en_x, coords_en_y, 'o:r')
plt.grid() # Esto me deja poner una grilla detrás!
plt.show()
Vemos que los puntos que se grafican son los $(x,y)$ = $\{(0,2);\,(2,4);\,(6,6);\,(8,14);\,(10,14)\}$
Dar formato y estética a los gráficos¶
En los gráficos anteriores cuando usamos plt.plot agregamos un término opcional para cambiar el color y la forma en la que se muestran los puntos que graficamos. Si usamos help(plt.plot) vemos que este argumento opcional es el [fmt] y que toma distintos strings con los que se puede cambiar la forma en que se ve el gráfico. La estructura general es
fmt = '[color][marker][line]'
donde marker hace referencia a la forma de los puntitos y line hace referencia a la forma de la linea.
Los strings válidos son un montón y los pueden ver en la documentación usando el help o en la documentación online de la función.
Acá dejamos algunas de las opciones para que se den una idea:
- color: 'r' (rojo), 'b' (azul), 'g' (verde), 'm' (magenta), etc
- marker: '.' (puntito chiquito), 'o' (punto grande), '^' (triángulo), '*' (estrella), etc
- line: '-' (sólida), '--' (de a rayas), ':' (de a puntos)
Título, labels y grid¶
Veamos también algunas funciones que nos permiten agregar información a la figura sobre la que estamos trabajando. Veamos un ejemplo con la función de los Arctic Monkeys.
import numpy as np
import matplotlib.pyplot as plt
# Definimos la función
def f(x):
return np.sin(x)*np.sin(20*x)
# Definimos el dominio y la imagen
x = np.linspace(0, 2*np.pi, 500)
y = f(x)
# Graficamos la función con linea sólida y le ponemos una etiqueta
plt.plot(x, y, '-', label='Arctic Monkeys')
# Agregamos título y etiquetamos los ejes
plt.title('Logo de los Arctic Monkeys', fontsize=16)
plt.xlabel('Eje de $\hat x$')
plt.ylabel('Eje de $\hat y$')
plt.grid()
plt.legend(loc="upper right")
plt.show()
La mayoría de las funciones son bastante descriptivas. Veamos solo algunas aclaraciones:
- El fontsize en el
plt.titlees un argumento opcional y también se puede utilizar en elplt.xlabeloplt.ylabelpor ejemplo. - El
plt.grid()agrega una grid si no exite y la elimina si existe. Si se quiere forzar a que aparezca se puede utilizarplt.grid(True) - El
plt.legend()es el que permite que se muestren los labels asociados a las lineas (en este caso es el "Arctic Monkeys"). Esta función tiene algunos parámetros opcionales que pueden probar como el fontsize o el loc. (verhelp(plt.legend)) - Podemos usar $\LaTeX$ en los labels!
Aca abajo dejamos una figura con los nombres de las distintas partes de una figura para que tengan a mano a la hora de googlear como cambiar cada parte. (fuente: este libro hermoso)
Ejercicio 2¶
En la misma figura:
- Graficar los datos cargados en las variables
t_teoricoey_teorico. Etiquetar los datos como 'Teorico'. - Graficar los datos cargados en las variables
t_laboey_labo. Etiquetar los datos como 'Laboratorio'. - Agregar la leyenda, etiquetar los ejes y poner
un título (pueden usar $\LaTeX$ si conocen la syntaxis)
Hacer ajustes¶
Vamos a levantar nuevamente unos datos para graficarlos igual que lo hacíamos antes.
import numpy as np
import matplotlib.pyplot as plt
# Importamos los datos
data = np.loadtxt('data_exp.csv', delimiter=',', skiprows=1, unpack=True)
# Los ponemos en variables
tiempo = data[0]
num_conejos = data[1]
err_num_conejos = data[2]
# O lo que es lo mismo
tiempo, num_conejos, err_num_conejos = data
# Graficamos los datos crudos
plt.plot(tiempo, num_conejos, 'ok', label='Datos')
# Labels, grilla, leyenda y show
plt.xlabel("Tiempo")
plt.ylabel("Número de conejos")
plt.grid()
plt.legend()
plt.show()
Los datos parecen tener un comportamiento exponencial. Queremos ver que tan exponencial es su distribución realizando un ajuste. Realizamos un ajuste primero definiendo la función que queremos comparar a los datos, $f$, con los parametros libres que nos interesa determinar, $$ \text{variable dependiente} = f(\text{variable independiente}, \text{parametros a determinar}). $$
En este caso, vamos a intentar con: $$ f(t, A, C) = Ae^{Ct} $$
Luego, para realizar un ajuste podemos usar la función curve_fit() que nos ofrece la librería/sublibrería scipy.optimize. Veamos cómo funciona con un ejemplo:
# Importo la función desde la librería
from scipy.optimize import curve_fit
# Declaro la función con la que quiero ajustar con parámetros genéricos
def f_ajuste(t, A, C):
exponencial = A*np.exp(C*t)
return exponencial
# Utilizo curve_fit() para el ajuste
popt, pcov = curve_fit(f_ajuste, tiempo, num_conejos, sigma = err_num_conejos)
# Imprimo en pantalla los valores de popt y pcov
print(f'Parámetros óptimos para A y C (popt): {popt}')
print(f'\nMatriz de covariancia de popt (pcov):\n{pcov}')
Parámetros óptimos para A y C (popt): [0.97349007 0.10119234] Matriz de covariancia de popt (pcov): [[ 4.92276135e-03 -1.89187395e-04] [-1.89187395e-04 7.51566964e-06]]
# Graficamos los DATOS con el error asociado
plt.errorbar(tiempo, num_conejos, yerr=err_num_conejos, fmt='ok', label='Datos',
capsize=2, capthick=1)
# Graficamos el AJUSTE
A, C = popt
# Variables finas para evaluar la funcion en más puntos que los medidos
x_ajuste = np.linspace(min(tiempo), max(tiempo), 1000)
y_ajuste = f_ajuste(x_ajuste, A, C)
plt.plot(x_ajuste, y_ajuste, '-r', label='Ajuste')
# Labels, grilla, leyenda y show
plt.xlabel("Tiempo")
plt.ylabel("Número de conejos")
plt.grid()
plt.legend()
plt.show()
Vemos que arriba la función recibe cinco inputs (3 obligatorios y 2 opcionales):¶
f: El primer argumento es la función modelo con la que queremos ajustar.La función modelo tiene que tener como primer parámetro la variable independiente (que corresponde por ejemplo al eje $\hat x$). Luego se usan los otros parámetros para poner las constantes que queremos hallar (en nuestro caso
AyC).xdata, ydata: Los otros dos argumentos son nuestros datosEstos son los datos crudos, los que usamos para visualizarlos normalmente.
sigma: El error asociado a cada valor enydataEs un parámetro opcional, para que el ajuste considere el error en la variable $y$.
Outputs de curve_fit()¶
popt: Parámetros Óptimos (array 1Dde numpy)Vemos que
curve_fit()nos devuelve como primer parámetro una lista con los valores óptimos para nuestro ajuste. En este caso son dos porque en la función modelo pusimos dos constantes (AyC). Estos valores vienen ordenados del mismo modo que están en la función modelo por lo quepopt[0]es elAóptimo ypopt[1]es elCóptimo.pcov: Covarianza de los Parámetros (array 2Dde numpy)Esta es la matriz de covarianza de los resultados que nos da
popt. Esta matriz nos da una idea de cuál es el error de este ajuste y de cuán ligados están estos errores entre sí. No es el objetivo del curso dar una clase de estadísitica, pero si conocemos que las variables son independientes podemos obtener el error de nuestros parámetros de ajuste como la raíz cuadrada de la covarianza. Para este ejemplo seríaerr_A = np.sqrt(pcov[0,0]) err_C = np.sqrt(pcov[1,1])
o
err_A, err_C = np.sqrt(np.diag(pcov))
Una cosita más: En este caso pusimos var_x_ajuste = var_x al momento de graficar, lo que implica que nuestro grafico del ajuste tiene 20 puntos (los mismos que var_x), pero como conocemos la función podríamos hacer un linspace del estilo np.linspace(min(var_x), max(var_x), 1000) para tener una mejor densidad de puntos. Esto esta bueno para hacer que los gráficos se vean menos acartonados, pero SOLAMENTE se puede usar a la hora de graficar.
Barras de error¶
- Si queremos poner barras de error a nuestro gráfico podemos usar
plt.errorbar(), que funciona de forma muy similar aplt.plot()pero nos permite dar un array con los errores o un error constante para todos.
Guardar un gráfico¶
Una vez tenemos un gráfico podemos guardarlo en el mismo directorio donde tenemos el archivo con el que estamos trabajando utilizando la linea
plt.savefig('nombre_del_archivo.png')
Ejercicio 3¶
Con los datos levantados en las variables t_labo e y_labo, realizar el ajuste de la función:
Donde $a = y_0$ y $b = -\frac{g}{2}$.
Agregar al gráfico del Ejercicio 3 este nuevo ajuste, junto con su label 'Ajuste'.
Guardar la imagen, e imprimir la gravedad en pantalla con su error.
Acá termina nuestra clase.¶
Si aún hay tiempo abajo tenemos armado algunos párrafos sobre distintas cosas piolas que podemos hacer. La elección esta en ustedes. (Si, se convirtió en un elija su propia aventura, no se la esperaba nadie esa)
- Pandas y DataFrames
- Histogramas
- Sistemas de ecuaciones diferenciales
- Algebra Lineal
- Calculo Simbólico
- Integración Numérica
- Funciones Locas
Pandas y DataFrames¶
Pandas es una librería bastante utilizada para Análisis de Datos y tiene demasiadas cosas disponibles, demasiadas. Aca en el curso no vamos a cubrir casi nada, sólo veremos algunas cosas.
Pandas, tal como NumPy y SymPy, trae un nuevo tipo de dato muy genial, los DataFrames. Recomendado es el tutorial de 10 minutos que aporta la misma página de Pandas, en los que hace un paneo sobre las cosas básicas de la librería. Hay un ejemplo práctico hecho por la FIFA en el que se trabaja un poco los datos de inscripción de una jornada anterior del cursito este.
Los DataFrames, el nuevo tipo de dato introducido, son básicamente tablas de datos donde cada columna tiene un nombre descriptivo. Podemos operar con las columnas por su nombre, y también elegir filas como haciamos con arrays antes.
Vamos a trabajar con un csv que pueden descargar de nuestro repositorio:
# Importamos pandas con un sobrenombre, como haciamos con numpy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# Para cargar un csv usamos read_csv() de pandas.
# El argumento es el nombre del archivo
df = pd.read_csv('datos_jugadores.csv')
# df es un "DataFrame", como si fuera una tabla de Excel cargada en Python.
display(df) # Display es como un "print" especial para cosas que no son texto.
| Equipo | Puesto | Jugador | Edad | Nacimiento | dia_mes | nro_mes | año | Altura_cm | Ciudad_nac | Pais_nac | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Tigre | MED | Aaron Molinas | 22 | 8/2/2000 | 02-ago | 8 | 2000 | 175.0 | NaN | Argentina |
| 1 | Banfield | DEF | Aarón Quirós | 21 | 10/31/2001 | 31-oct | 10 | 2001 | NaN | NaN | Argentina |
| 2 | Tigre | DEF | Abel Luciatti | 30 | 2/18/1993 | 18-feb | 2 | 1993 | 178.0 | Morón | Argentina |
| 3 | Velez | DEL | Abiel Osorio | 20 | 6/13/2002 | 13-jun | 6 | 2002 | 182.0 | NaN | Argentina |
| 4 | San Lorenzo | DEL | Adam Bareiro | 26 | 7/26/1996 | 26-jul | 7 | 1996 | 184.0 | Asuncion | Paraguay |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 639 | Atl Tucuman | DEF | Wilson Ibarrola | 26 | 7/2/1996 | 02-jul | 7 | 1996 | 172.0 | NaN | Paraguay |
| 640 | Sarmiento (J) | MED | Yair Arismendi | 25 | 4/5/1998 | 05-abr | 4 | 1998 | NaN | NaN | Argentina |
| 641 | Union | MED | Yeison Gordillo | 30 | 6/25/1992 | 25-jun | 6 | 1992 | 176.0 | Miranda (Cauca) | Colombia |
| 642 | Atl Tucuman | DEF | Yonathan Cabral | 30 | 5/10/1992 | 10-may | 5 | 1992 | 188.0 | Isidro Casanova | Argentina |
| 643 | Estudiantes (LP) | DEF | Zaid Romero | 23 | 12/15/1999 | 15-dic | 12 | 1999 | 193.0 | Mendoza | Argentina |
644 rows × 11 columns
# Obtener cantidad de filas y columnas total:
print(df.shape)
# Obtener nombres de las columnas:
print(df.columns)
(644, 11)
Index(['Equipo', 'Puesto', 'Jugador', 'Edad', 'Nacimiento', 'dia_mes',
'nro_mes', 'año', 'Altura_cm', 'Ciudad_nac', 'Pais_nac'],
dtype='object')
Nos gustaría operar de alguna manera con esta tabla, para hacer algun tipo de análisis de los datos.
Preguntas posibles:
- Cuál es la distribución de edad/altura/año de nacimiento?
- Cuantos jugadores de cada País/Puesto/mes de nacimiento hay?
- Hay dos jugadores que compartan cumpleaños en un mismo equipo?
Cuál es la distribución de edad?¶
Para esta pregunta vamos a armar un histograma de los valores de edad.
# Para seleccionar una columna usamos:
altura = df['Altura_cm']
altura.hist()
<Axes: >
# Dejemos un poco más lindo el gráfico
altura.hist(rwidth=0.9, edgecolor='black', color='lightblue')
# Label de los ejes
plt.xlabel('Altura [cm]', fontsize=13)
plt.ylabel('Ocurrencias', fontsize=13)
# Apagamos la grid que viene por defecto y despues colocamos la nuestra
plt.grid()
plt.grid(ls='--', axis='y', color='gray')
plt.title('Distribución de altura de los jugadores')
plt.show()
Como hacemos si queremos crear una columna nueva? Por ejemplo, el logaritmo de la altura?
df['Altura_log'] = np.log(df['Altura_cm'])
# Dejemos un poco más lindo el gráfico
df['Altura_log'].hist(rwidth=0.9, edgecolor='black', color='lightblue')
# Label de los ejes
plt.xlabel('log(Altura)', fontsize=13)
plt.ylabel('Ocurrencias', fontsize=13)
# Apagamos la grid que viene por defecto y despues colocamos la nuestra
plt.grid()
plt.grid(ls='--', axis='y', color='gray')
plt.title('Distribución de altura de los jugadores')
plt.show()
Cuántos jugadores de cada pais hay?¶
# Podemos usar la función groupby de pandas para agrupar todas las filas
# con el mismo país de nacimiento.
df_gb = df.groupby('Pais_nac')
display(df_gb)
<pandas.core.groupby.generic.DataFrameGroupBy object at 0x7f43887eb430>
Vemos que esto asi solo no nos dice absolutamente nada. Una vez que tenemos agrupado por una columna nuestro DataFrame, podemos hacer cálculos como el promedio, la suma de los valores, la cantidad de veces que aparece una fila con el valor que agrupamos, el máximo, etc. A estas se las conoce como funciones de agregación.
# Para contar la cantidad de filas seleccionamos una columna cualquiera y usamos 'count'.
df_gb['Equipo'].count()
Pais_nac Argentina 544 Chile 5 Colombia 22 Ecuador 3 Espana 3 Mexico 1 Paraguay 28 Peru 2 Portugal 1 Uruguay 33 Venezuela 2 Name: Equipo, dtype: int64
# Si quisieramos el máximo de altura, por ejemplo:
df_gb['Altura_cm'].max()
Pais_nac Argentina 194.0 Chile 180.0 Colombia 190.0 Ecuador 176.0 Espana 192.0 Mexico 178.0 Paraguay 188.0 Peru 185.0 Portugal 176.0 Uruguay 190.0 Venezuela 186.0 Name: Altura_cm, dtype: float64
Hay dos jugadores que compartan cumpleaños en un mismo equipo?¶
# Para responder esto, podemos contar cuantas veces se repite
# el mismo valor de "Equipo" y "dia_mes". Para hacer esto, podemos aprovechar
# a usar devuelta el groupby, y contar cuantas veces está esa combinaciones
# de valores
df_gb = df.groupby(['Equipo', 'dia_mes'])['Puesto'].count()
display(df_gb)
Equipo dia_mes
Argentinos 01-dic 1
02-feb 1
04-jul 1
07-nov 1
08-ago 1
..
Velez 25-feb 1
26-jun 1
28-abr 1
30-ene 1
31-ene 1
Name: Puesto, Length: 622, dtype: int64
# Ahora podemos filtrar usando esta sintaxis
df_filtered = df_gb[df_gb > 1]
display(df_filtered)
Equipo dia_mes
Argentinos 25-ene 2
Atl Tucuman 10-may 2
17-feb 2
Banfield 01-feb 2
Belgrano 29-ago 2
Boca Juniors 04-jul 2
22-feb 2
23-jul 2
24-abr 2
Central Cba (SdE) 13-ene 2
27-oct 2
Colon 05-sep 2
13-ene 2
Def y Justicia 08-feb 2
Estudiantes (LP) 29-mar 2
Gimnasia (LP) 05-may 2
Godoy Cruz 02-mar 2
16-feb 2
Lanus 22-mar 2
River Plate 29-mar 2
Tigre 07-feb 2
Union 30-jul 2
Name: Puesto, dtype: int64
# Como hacemos un filtro un poco mas complicado?
# Elijamos las personas que nacieron en Agosto, en Argentina y su altura es mayor que 150 cm
df_filtered = df[(df['nro_mes'] == 8) & (df['Altura_cm'] > 150) & (df['Pais_nac'] == 'Argentina')]
display(df_filtered)
| Equipo | Puesto | Jugador | Edad | Nacimiento | dia_mes | nro_mes | año | Altura_cm | Ciudad_nac | Pais_nac | Altura_log | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Tigre | MED | Aaron Molinas | 22 | 8/2/2000 | 02-ago | 8 | 2000 | 175.0 | NaN | Argentina | 5.164786 |
| 67 | Independiente | DEL | Braian Martínez | 23 | 8/18/1999 | 18-ago | 8 | 1999 | 170.0 | NaN | Argentina | 5.135798 |
| 146 | Belgrano | DEF | Erik Godoy | 29 | 8/16/1993 | 16-ago | 8 | 1993 | 185.0 | Buenos Aires | Argentina | 5.220356 |
| 202 | Platense | MED | Franco Diaz | 22 | 8/28/2000 | 28-ago | 8 | 2000 | 180.0 | NaN | Argentina | 5.192957 |
| 219 | Belgrano | MED | Gabriel Compagnucci | 31 | 8/29/1991 | 29-ago | 8 | 1991 | 178.0 | Monte Buey | Argentina | 5.181784 |
| 246 | Def y Justicia | MED | Gonzalo Castellani | 35 | 8/10/1987 | 10-ago | 8 | 1987 | 183.0 | Buenos Aires | Argentina | 5.209486 |
| 248 | Central Cba (SdE) | DEF | Gonzalo Goñi | 24 | 8/16/1998 | 16-ago | 8 | 1998 | 187.0 | Pedro Luro | Argentina | 5.231109 |
| 265 | Newells | DEF | Guillermo Ortíz | 30 | 8/9/1992 | 09-ago | 8 | 1992 | 183.0 | Rosario | Argentina | 5.209486 |
| 279 | Atl Tucuman | MED | Ignacio Maestro Puch | 19 | 8/13/2003 | 13-ago | 8 | 2003 | 180.0 | San Miguel de Tucuman | Argentina | 5.192957 |
| 282 | Platense | MED | Ignacio Schor | 22 | 8/4/2000 | 04-ago | 8 | 2000 | 183.0 | Buenos Aires | Argentina | 5.209486 |
| 343 | Rosario Central | DEF | Juan Cruz Komar | 26 | 8/13/1996 | 13-ago | 8 | 1996 | 190.0 | Rosario | Argentina | 5.247024 |
| 368 | Talleres (C) | DEF | Julio Buffarini | 34 | 8/18/1988 | 18-ago | 8 | 1988 | 172.0 | Gral. Cabrera | Argentina | 5.147494 |
| 410 | Huracan | ARQ | Lucas Chaves | 27 | 8/9/1995 | 09-ago | 8 | 1995 | 179.0 | Martín Coronado | Argentina | 5.187386 |
| 416 | Velez | DEL | Lucas Janson | 28 | 8/16/1994 | 16-ago | 8 | 1994 | 171.0 | NaN | Argentina | 5.141664 |
| 424 | Argentinos | DEF | Lucas Villalba | 28 | 8/19/1994 | 19-ago | 8 | 1994 | 177.0 | NaN | Argentina | 5.176150 |
| 442 | Belgrano | ARQ | Manuel Vicentini | 32 | 8/29/1990 | 29-ago | 8 | 1990 | 187.0 | Santa Fé | Argentina | 5.231109 |
| 459 | Tigre | DEF | Martín Ortega | 23 | 8/20/1999 | 20-ago | 8 | 1999 | 183.0 | NaN | Argentina | 5.209486 |
| 467 | Talleres (C) | DEF | Matías Catalán | 30 | 8/19/1992 | 19-ago | 8 | 1992 | 172.0 | Mar del Plata | Argentina | 5.147494 |
| 487 | Central Cba (SdE) | MED | Mauro Pittón | 28 | 8/8/1994 | 08-ago | 8 | 1994 | 174.0 | NaN | Argentina | 5.159055 |
| 502 | Argentinos | DEF | Miguel Torrén | 34 | 8/12/1988 | 12-ago | 8 | 1988 | 179.0 | Santa Fé | Argentina | 5.187386 |
| 504 | Banfield | DEL | Milton Gimenez | 26 | 8/12/1996 | 12-ago | 8 | 1996 | 184.0 | Grand Bourg | Argentina | 5.214936 |
| 545 | Newells | MED | Pablo Pérez | 37 | 8/10/1985 | 10-ago | 8 | 1985 | 179.0 | Rosario | Argentina | 5.187386 |
| 571 | Argentinos | DEL | Rodrigo Cabral | 22 | 8/8/2000 | 08-ago | 8 | 2000 | 170.0 | Mercedes | Argentina | 5.135798 |
| 590 | Instituto | DEL | Santiago Rodríguez | 25 | 8/23/1997 | 23-ago | 8 | 1997 | 171.0 | San Luis | Argentina | 5.141664 |
| 598 | Union | ARQ | Sebastián Moyano | 32 | 8/26/1990 | 26-ago | 8 | 1990 | 188.0 | Mendoza | Argentina | 5.236442 |
| 629 | Huracan | MED | Valentín Burgoa | 22 | 8/16/2000 | 16-ago | 8 | 2000 | 175.0 | Guaymallén | Argentina | 5.164786 |
| 635 | Velez | DEL | Walter Bou | 29 | 8/25/1993 | 25-ago | 8 | 1993 | 176.0 | Concordia | Argentina | 5.170484 |
# Si quiero acceder a la fecha de nacimiento, del primer jugador de los que tenemos filtrados, como puedo hacer?
# Primero el numero de fila, luego la columna a la que queremos entrar
df_filtered.loc[0, 'dia_mes'] = 'Esto es una prueba'
display(df_filtered)
| Equipo | Puesto | Jugador | Edad | Nacimiento | dia_mes | nro_mes | año | Altura_cm | Ciudad_nac | Pais_nac | Altura_log | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Tigre | MED | Aaron Molinas | 22 | 8/2/2000 | Esto es una prueba | 8 | 2000 | 175.0 | NaN | Argentina | 5.164786 |
| 67 | Independiente | DEL | Braian Martínez | 23 | 8/18/1999 | 18-ago | 8 | 1999 | 170.0 | NaN | Argentina | 5.135798 |
| 146 | Belgrano | DEF | Erik Godoy | 29 | 8/16/1993 | 16-ago | 8 | 1993 | 185.0 | Buenos Aires | Argentina | 5.220356 |
| 202 | Platense | MED | Franco Diaz | 22 | 8/28/2000 | 28-ago | 8 | 2000 | 180.0 | NaN | Argentina | 5.192957 |
| 219 | Belgrano | MED | Gabriel Compagnucci | 31 | 8/29/1991 | 29-ago | 8 | 1991 | 178.0 | Monte Buey | Argentina | 5.181784 |
| 246 | Def y Justicia | MED | Gonzalo Castellani | 35 | 8/10/1987 | 10-ago | 8 | 1987 | 183.0 | Buenos Aires | Argentina | 5.209486 |
| 248 | Central Cba (SdE) | DEF | Gonzalo Goñi | 24 | 8/16/1998 | 16-ago | 8 | 1998 | 187.0 | Pedro Luro | Argentina | 5.231109 |
| 265 | Newells | DEF | Guillermo Ortíz | 30 | 8/9/1992 | 09-ago | 8 | 1992 | 183.0 | Rosario | Argentina | 5.209486 |
| 279 | Atl Tucuman | MED | Ignacio Maestro Puch | 19 | 8/13/2003 | 13-ago | 8 | 2003 | 180.0 | San Miguel de Tucuman | Argentina | 5.192957 |
| 282 | Platense | MED | Ignacio Schor | 22 | 8/4/2000 | 04-ago | 8 | 2000 | 183.0 | Buenos Aires | Argentina | 5.209486 |
| 343 | Rosario Central | DEF | Juan Cruz Komar | 26 | 8/13/1996 | 13-ago | 8 | 1996 | 190.0 | Rosario | Argentina | 5.247024 |
| 368 | Talleres (C) | DEF | Julio Buffarini | 34 | 8/18/1988 | 18-ago | 8 | 1988 | 172.0 | Gral. Cabrera | Argentina | 5.147494 |
| 410 | Huracan | ARQ | Lucas Chaves | 27 | 8/9/1995 | 09-ago | 8 | 1995 | 179.0 | Martín Coronado | Argentina | 5.187386 |
| 416 | Velez | DEL | Lucas Janson | 28 | 8/16/1994 | 16-ago | 8 | 1994 | 171.0 | NaN | Argentina | 5.141664 |
| 424 | Argentinos | DEF | Lucas Villalba | 28 | 8/19/1994 | 19-ago | 8 | 1994 | 177.0 | NaN | Argentina | 5.176150 |
| 442 | Belgrano | ARQ | Manuel Vicentini | 32 | 8/29/1990 | 29-ago | 8 | 1990 | 187.0 | Santa Fé | Argentina | 5.231109 |
| 459 | Tigre | DEF | Martín Ortega | 23 | 8/20/1999 | 20-ago | 8 | 1999 | 183.0 | NaN | Argentina | 5.209486 |
| 467 | Talleres (C) | DEF | Matías Catalán | 30 | 8/19/1992 | 19-ago | 8 | 1992 | 172.0 | Mar del Plata | Argentina | 5.147494 |
| 487 | Central Cba (SdE) | MED | Mauro Pittón | 28 | 8/8/1994 | 08-ago | 8 | 1994 | 174.0 | NaN | Argentina | 5.159055 |
| 502 | Argentinos | DEF | Miguel Torrén | 34 | 8/12/1988 | 12-ago | 8 | 1988 | 179.0 | Santa Fé | Argentina | 5.187386 |
| 504 | Banfield | DEL | Milton Gimenez | 26 | 8/12/1996 | 12-ago | 8 | 1996 | 184.0 | Grand Bourg | Argentina | 5.214936 |
| 545 | Newells | MED | Pablo Pérez | 37 | 8/10/1985 | 10-ago | 8 | 1985 | 179.0 | Rosario | Argentina | 5.187386 |
| 571 | Argentinos | DEL | Rodrigo Cabral | 22 | 8/8/2000 | 08-ago | 8 | 2000 | 170.0 | Mercedes | Argentina | 5.135798 |
| 590 | Instituto | DEL | Santiago Rodríguez | 25 | 8/23/1997 | 23-ago | 8 | 1997 | 171.0 | San Luis | Argentina | 5.141664 |
| 598 | Union | ARQ | Sebastián Moyano | 32 | 8/26/1990 | 26-ago | 8 | 1990 | 188.0 | Mendoza | Argentina | 5.236442 |
| 629 | Huracan | MED | Valentín Burgoa | 22 | 8/16/2000 | 16-ago | 8 | 2000 | 175.0 | Guaymallén | Argentina | 5.164786 |
| 635 | Velez | DEL | Walter Bou | 29 | 8/25/1993 | 25-ago | 8 | 1993 | 176.0 | Concordia | Argentina | 5.170484 |
Lambda functions¶
En Python existe otra forma de declarar funciones además de la que vimos anteriormente, que es utilizando la palabra lambda. La función promediar_dos_numeros que vimos más temprano se escribiría de la siguiente forma
promedio = lambda num1, num2: (num1 + num2) / 2
Los lambda son funciones "anónimas", y fueron introducidas a Python para código que sigue un estilo de programación llamado "funcional". Generalmente, la idea no es asignarla a una variable, sino usarla directamente dentro de otra función. Por ejemplo, la función filter:
pares = list(filter(lambda x: x % 2 == 0, range(10)))
Lo introducimos por completitud, por si se la encuentran leyendo el código de otra persona, pero lo recomendado es definir funciones "normales" con def.
Histogramas¶
Una cosita más que nos va a ser útil a la hora de dejar el Oriyin sin instalar es poder hacer histogramas. Con pyplot eso lo podemos obtener de la función hist.
Recordemos que en un histograma 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.
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 sintéticos usando random.normal de NumPy . Esto de crear datos lo hacemos acá a modo de ejemplo, en la vida real uno importaria algun dataset de las formas que ya hemos visto.
import numpy as np
import matplotlib.pyplot as plt
# Distribución gaussiana centrada en 100, con 15 de desviación estándar
data = np.random.normal(100, 15, 2000)
n, bins, patches = plt.hist(data, bins=20, edgecolor='black', facecolor='royalblue', alpha=0.75)
# n 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
# OBS1: Si no saben qué cantidad de bines elegir, pueden dejarle la
# eleccion a matplotlib, usando: bins='auto'
# OBS2: Si quieren forzar bines específicos, pueden pasarle un array o lista
# con los comienzos de cada bin. Ej:
#
# Esto da bines de ancho 0.1:
# bins = [1, 1.1, 1.2, 1.3, 1.4, 1.5]
#
# Indistintamente pueden usar cosas que ya vimos para generar ararys:
# bins = np.arange(1, 1.6, 0.1)
print(f'n:\n{n}\n')
print(f'bins:\n{bins}')
n: [ 2. 0. 2. 12. 24. 55. 101. 190. 245. 279. 291. 291. 220. 137. 77. 46. 18. 6. 1. 3.] bins: [ 42.26532864 47.92479703 53.58426543 59.24373383 64.90320223 70.56267062 76.22213902 81.88160742 87.54107581 93.20054421 98.86001261 104.51948101 110.1789494 115.8384178 121.4978862 127.1573546 132.81682299 138.47629139 144.13575979 149.79522819 155.45469658]
# Si lo quisieramos normalizar el area hay que agregar una opcion mas que es density=True como para
# que entienda que queremos ver la "densidad de probabilidad", forma cheta para decir: normalizar el area.
n, bins, patches = plt.hist(data, bins=20, edgecolor='black', facecolor='royalblue', alpha=0.75, density=True)
Y ya que estamos, para concientizar acerca de los peligros a la hora de la elección de bins, graficamos algunos histogramas superpuestos.
n, bins, patches = plt.hist(data, bins=100, density=True, edgecolor='black', facecolor='royalblue', alpha=0.75)
n, bins, patches = plt.hist(data, bins=20 , density=True, edgecolor='black', facecolor='purple' , alpha=0.5)
n, bins, patches = plt.hist(data, bins=2 , density=True, edgecolor='black', facecolor='indianred' , alpha=0.3)
Ejercicio 6¶
La función randn que usamos nos brinda números aleatorios distribuidos de manera gaussiana (distribución normal).
Con lo mostrado arriba, cree un histograma con una distribución normal (igual que como hicimos ya) y superpongale una curva teórica de una gausiana.
AYUDA: Busquen "plot normal distribution python", en particular busquen la librería
scipy.Haga el mismo ejercicio de recién pero con números aleatorios distribuidos de manera uniforme
Lo mismo pero con números con una distribución pareto, o la distribución que prefieran.
Info: https://docs.scipy.org/doc/numpy-1.14.0/reference/routines.random.html
Resolución de un sistema de ecuaciones diferenciales ordinarias (ODE)¶
Vamos con un caso simple y conocido por la mayoría: el infaltable problema del péndulo. Arrancamos desde donde sabemos todos. (Donde sabemos todos? Si no saben de donde salió esto no se hagan problema, es solo algo físico, no lo podemos evitar)
$$ \frac{d^2\theta}{dt^2} + \omega^2 \theta = 0$$con $\omega^2 = \frac{g}{l}$
Para resolver numéricamente este problema, proponemos utilizar la función odeint de la biblioteca scipy.integrate (Otra opción podría ser la función solve_ivp, para problemas de valores iniciales), y esa función utiliza un método no taaaaan distinto al método de Euler para la solución de ODE's, (inserte conocimientos de Cálculo Numérico aquí, sí, dale, cursala antes de recibirte) por lo que, para un problema con derivadas de segundo orden debemos armar un sistema de ecuaciones de primer orden.
Sea $\phi = \dot{\theta} \rightarrow \dot{\phi} = \ddot{\theta}$
Nos queda entonces
o bien
o también
$$ \dot{\vec{X}} = A \vec{X} $$Escencialmente, nuestro odeint va a intentar resolver ese sistema para distintos $t$ mientras le hayamos dado un valor inicial de donde comenzar. Nada demasiado extraño. Así que lo que la función va a necesitar es una función que tome como primer argumento $\vec{X}$, segundo argumento $t$ (por ser la variable independiente) y luego los demás parametros que necesite (en este caso sean $g$ y $l$) y opere para obtener $\dot{\vec{X}}$.
from scipy.integrate import odeint
import numpy as np
def ecdif(X, t, g, l):
theta, phi = X
omega2 = g/l
return [phi, -omega2 * theta]
g = 9.8
l = 2
X0 = [np.pi/4, 0] # Inicio a 45 grados con velocidad = 0
t = np.linspace(0, 10, 101)
solucion = odeint(ecdif, X0, t, args = (g,l))
plt.plot(t,solucion[:,0], label = 'theta')
plt.plot(t,solucion[:,1], label = 'phi')
plt.title('Resolucion del pendulo')
plt.xlabel('Tiempo')
plt.ylabel('Valores')
plt.grid(True)
plt.legend(loc = 'best')
plt.show()
Si quieren ver otro ejemplo que esto, recomendamos el resuelto que hizo un Ex-FIFA del oscilador con forzante sin aproximación y comparado con la solución analítica (sí, esa que sin aproximación no buscó nadie). También recomendamos, a cuento de esto, un pasito más en resolución de ODE's que es una simulación de uno de los ejercicios de F1, donde se arma una animación y una visualización más interactiva con la biblioteca ipywidgets de tantas otras posibles.
Ejercicio¶
- Resuelvan el oscilador armónico amortiguado con el parámetro $\gamma$
- Vuelvan a todas esas guías que nadie terminó, busquen los ejercicios con asterisco, métanlos en Python y aprecien el poder del cálculo numérico.
Un poco de algebra lineal con NumPy¶
NumPy trae muchas funciones para resolver problemas típicos de algebra lineal usando a los arrays como vectores. El que nos interesa en general es el de autovalores y autovectores y el sistema de ecuaciones lineales, pero empecemos con un ejemplo más fácil:
# cargamos el módulo de algebra lineal
from numpy import linalg
v = np.array([1, 1, 1])
w = np.array([2, 2, 2])
z = np.array([1, 0, 1])
norma =linalg.norm(v) # la norma 2 / módulo del vector v
print(norma)
print(np.sqrt(v[0]**2 + v[1]**2 + v[2]**2)) # calculado a mano == sqrt(3)
print(norma == np.sqrt(3)) # y numpy sabe que son lo mismo
1.7320508075688772 1.7320508075688772 True
Ahora si, usemos los vectores que creamos recién para crear una matriz y digamosle a NumPy que calcule los autovectoresy autovalores de esa matriz:
matriz = np.array([v, w, z], dtype=np.float64)
#eig devuelve una tupla de arrays con los autovalores en un array 1D y los autovec en un array 2D
eigens = linalg.eig(matriz)
autvals, autvecs = eigens
print('Los autovalores:', autvals)
print()
print('Los autovectores:', autvecs)
#se terminó el problema
Los autovalores: [ 3.41421356e+00 -1.23150120e-16 5.85786438e-01] Los autovectores: [[ 4.39732612e-01 7.07106781e-01 -3.03890631e-01] [ 8.79465224e-01 -1.09165123e-16 -6.07781262e-01] [ 1.82143212e-01 -7.07106781e-01 7.33656883e-01]]
Y para un sistema de ecuaciones del tipo $Ax = b$:
mat = np.array([[1, 2, 5], [2, 5, 8], [4, 0, 8]], dtype=np.float64)
b = np.array([1, 2, 3])
x = linalg.solve(mat, b) # resuelve el sistema A*x = b
print(x)
# se terminó el problema
[0.67857143 0.07142857 0.03571429]
Por supuesto, también se puede hacer producto matriz con vector, y... oh si, se pueden calcular inversas.
print(np.dot(mat,x))
print(linalg.inv(mat))
[1. 2. 3.] [[-1.42857143 0.57142857 0.32142857] [-0.57142857 0.42857143 -0.07142857] [ 0.71428571 -0.28571429 -0.03571429]]
Ejercicio¶
Parecen incrédulos. Fabriquen entonces la matriz $$ \begin{equation} A = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \end{equation} $$
cuyos autovalores son $\lambda_1 = 1$ y $\lambda_2 = -1 $, hallen sus autovalores y autovectores (a mano y con Python) y calculen $Ax$ con $$ \begin{equation} x = \begin{pmatrix} 1\\0 \end{pmatrix} \end{equation} $$
Ejercicio¶
Inventen una matriz de 5x5 (con el método que quieran y ¡que no sea la identidad!) y averigüen si es invertible. Si están prestando atención, saben que para resolver este problema les conviene ver la documentación de numpy.linalg.
Ayuda: ¿Es necesario calcular la inversa? ¿Se acuerdan algo de álgebra lineal?
Cálculo simbólico¶
Ahora no sólo nos vamos a emancipar del Oriyin, sino también del bendito Wolfram Alpha (Disclaimer: si bien la librería que se usa es muy buena, todavía no llega a tener el nivel de software como Mathematica, Maple, etc). Así es, Python nos va a permitir hacer cálculos simbólicos, como las integrales que nunca supimos calcular. Todo eso y más, en el paquete SymPy:
Hemos dicho que en general no es una buena práctica, pero este es uno de los pocos casos para los cuales se justifica importar toda la librería.
from sympy import *
from sympy.abc import *
Les daremos ahora un breve tour por algunas de las funciones que nos ofrece SymPy. Así como numpy nos introdujo el array como tipo de variable, las expresiones algebraicas en sympy son del tipo symbols. Estas pueden representar números enteros, reales, funciones, etc.
A diferencia de las librerías anteriores, en las cuales usualmente escribimos los comandos dentro de un archivo a ejectutar todo junto (script); a modo de demostración utilizaremos SymPy de forma interactiva. Esta es posiblemente la manera en la que usaron Wolfram Alpha o Mathematica (para quienes lo hayan hecho).
Vimos aquí que una cantidad de variables se definieron como symbols de distintos tipos. Veamos ahora algunas de las posibilidades que tenemos:
expand( (x + y)**2 )
factor( x**6 - 1 )
Utilizando variables simbólicas de tipo integer, podemos hacer algunas sumatorias, por ejemplo, la famosa $$\sum_{k=0}^{m}k$$
Sum(k, (k,0,m) ).doit().factor() #el comando .doit() evalúa la suma
O incluso algunas series, como $$ \sum_{n=1}^{\infty}\frac{1}{n^2} $$
Sum(1/n**2, (n, 1, oo)).doit() # el infinito se escribe como oo (dos o minúscula)
Posiblemente si queremos algo de cálculo simbólico, sea para calcular derivadas e integrales que nos molesten. Veamos algunos ejemplos $\dfrac{\text{d}x^n}{\text{d}x}$
diff(x**n , x).simplify()
diff(sin(tan(8*x)), x)
out = diff(sin(tan(8*x)), x).simplify()
# Podemos hacer sustituciones de simbolos por otros simbolos!!
out.subs(x, sqrt(y))
# Hasta podemos sustituir por valores numéricos y evaluar la funcion!!
print(f"{out.subs(x, 8)} = {out.subs(x, 8).evalf()}")
8*cos(tan(64))/cos(64)**2 = -36.5316624851842
Si alguno desconfia, puede verlo desde WolframAlfa
Otro ejemplo! Derivadas parciales: $\dfrac{\partial^2}{\partial y \partial x}\left( x\sin{y} \right)$
diff(x*sin(y), x, y) # qué pasa si en lugar de x,y ponemos y,y?
También se pueden hacer integrales indefinidas o definidas con dominio infinito/acotado, obviamente luego se podría reemplazar en la expresión final por algun valor numérico, si así se quisiera $$\int \dfrac{1}{1+x^2}\text{d}x$$
integrate(1/(1+x**2), x)
integrate( sin(x)/x, (x,-oo,oo) )
Algo importantisimo y fabuloso de Sympy (algo que sorprendentemente NO tiene Mathematica ni Matlab) es que a cualquier expresión, podemos envolverla en un latex() para que devuelva la misma expresión que renderizó, pero en el formato de $\LaTeX$. Derecho para mandar al informe o trabajo sin tipearla.
out = diff(sin(tan(8*x))*log(x), x)
out
print(f"Normal: {out}")
print(f"LaTeX: {latex(out)}")
Normal: (8*tan(8*x)**2 + 8)*log(x)*cos(tan(8*x)) + sin(tan(8*x))/x
LaTeX: \left(8 \tan^{2}{\left(8 x \right)} + 8\right) \log{\left(x \right)} \cos{\left(\tan{\left(8 x \right)} \right)} + \frac{\sin{\left(\tan{\left(8 x \right)} \right)}}{x}
Ejercicio 9¶
Calculen la siguiente integral: $$ \int_{-\infty}^{+\infty} e^{-x^2} \text{d}x $$ (opcional: verificar el resultado a mano!)
Integración numérica¶
Siempre es útil tener a mano una rutina de integración numérica. Puede ser para el caso en el cual la integral analítica sea muy complicada o no estándar. También es fundamental para poder integrar datos obtenidos experimentalmente, sin asumir alguna función que los modele. Para ambos casos, el paquete relevante será scipy.integrate
A partir de una función predefinida¶
En este caso, la función a importar es quad, porque usa un método de cuadraturas
from scipy.integrate import quad
Definamos una función para ser integrada. A modo de ejemplo, calcularemos
$$\int_a^b \dfrac{\text{e}^{-x^2 / 2}}{\sqrt{2\pi}}\text{d}x$$def integrando(x): return np.exp(-x**2 / 2)/np.sqrt(2*np.pi)
Ahora, llamamos a quad para integrar esta función entre $a=-1$ y $b=1$
a = -1
b = 1
integral, error = quad(integrando, a, b)
print(integral, error)
0.682689492137086 7.579375928402476e-15
Podemos, con este método, incluso obtener la primitiva numéricamente. Por ejemplo
$$F(x) = \int_a^x \dfrac{\text{e}^{-t^2 / 2}}{\sqrt{2\pi}}\text{d}t$$# Definimos este dominio
x = np.linspace(-5,5,100)
# Ahora obtenemos la primitiva
prim = np.zeros(x.size)
for i in range(x.size):
prim[i], error = quad(integrando, a, x[i])
Veamos la gráfica para este dominio
plt.figure()
plt.plot(x, integrando(x), 'b', label = 'integrando')
plt.plot(x, prim, 'r', label = 'primitiva')
plt.grid(True)
plt.legend()
<matplotlib.legend.Legend at 0x7f4386bb67a0>
Funciones especiales¶
Tal vez en algún momento de la vida nos encontremos con funciones más exóticas que el seno, coseno y exponenciales, funciones que las conocemos de nombre pero de graficarlas ni hablemos. Por suerte, la biblioteca scipy cuenta con abanico muy grande de estas funciones con nombre propio. En la documentación (bloque de ayuda del spyder) podemos encontrar la lista completa de funciones, buscando scipy.special
import scipy.special as sp
Empecemos graficando un ejemplo común en probabilidad, la función error:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
x = np.linspace(-5,5,100) # generamos el dominio
y = sp.erf(x) # llama a la función error de la librería
plt.plot(x,y, label='función error')
plt.grid(True)
plt.legend(loc='best')
plt.show()
Dentro de las muchas funciones que nos ofrece la librería, otro ejemplo importante en la física son las funciones de Bessel $J_{\nu}(x)$. Hay una función para cada valor del índice $\nu$ (que llamamos orden). En la librería, se llaman con el comando jv(i,x), donde la primer entrada corresponde al orden $\nu$.
Como demostración, veamos ahora algunos de sus gráficos:
num = 500
r = np.linspace(0,10,num)
for i in range(5):
y_i = sp.jv(i,r) # para cada i, grafica la función de orden i
plt.plot(r,y_i,label='orden {}'.format(i))
plt.legend(loc='best')
plt.grid(True)
plt.xlim(0, 10)
plt.title('Funciones de Bessel')
Text(0.5, 1.0, 'Funciones de Bessel')
Ejercicio 8¶
- Busquen en la documentación la función gamma, y grafiquenla en el intervalo $\left[-5,5\right]$
- (Opcional) Googleen a ver qué es lo tan importante de esa función.