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).
In [ ]:
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"
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()
In [ ]:
math.cos(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 (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:
In [3]:
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 = ', np.e)
#print(f"El numero e = {np.e}")
#print("El numero e = {}".format(np.e))
print('O el numero Pi = ', np.pi)
#print(f"O el numero Pi = {np.pi}")
#print("O el numero Pi = {}".format(np.pi))
# Podemos calcular senos y cosenos de numeros igual que con math!
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, 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).
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 [4]:
array_1 = np.array([1, 2, 3, 4]) # array toma como argumento un vector_like (lista, tupla, otro array...vector_like)
array_2 = np.array([5, 6, 7, 8])
# Se acuerdan como antes pasabamos de un objeto tipo 'range' a un objeto tipo 'list' porque eran "parecidos"?
# Aca hacemos lo mismo. Pasamos de un objeto tipo 'list' a un objeto tipo 'ndarray'
print("ARRAYS:")
print(f"{array_1} + {array_2} = {array_1 + array_2}") # aca los sumo
print(f"{array_1} * {array_2} = {array_1*array_2}") # acá los multiplico
# 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(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 con su índice:
In [5]:
print(array_1[0], array_1[1], array_1[2], array_1[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(array_2[-1]) # agarro al último elemento de b
print(array_2[0:3]) # desde el primero hasta el 3 (no incluido el final, nunca se incluye)
También podemos iterar sobre ellos con for , de las mismas formas que habíamos visto para iterar listas.
In [6]:
print("array_1:")
for num in array_1: # Aca itero sobre los ELEMENTOS de array_1
print(num)
print('\n') # Esto es solo para que agregue una linea nueva y no queden pegados
print("array_2:")
for j in range(len(array_2)): # Aca itero sobre los INDICES de array_2
print(array_2[j])
Para crear arrays de dos dimensiones (o más), podemos aplicar la función np.array sobre una lista de listas:
In [7]:
c = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) # cada lista corresponde a una fila de la matriz
print(c)
# Aca hacer c[0] dara la primer FILA de la matrix.
# Si queremos sacar la primer componente habra que hacer c[0, 0]
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 [8]:
# Equiespaciados >> Elijo CANTIDAD DE PUNTOS en un dado intervalo [a, b]
equi_linspace = np.linspace(0, 1, 10) # 9 números equiespaciados linealmente entre 0 y 1 (LINear SPACE = linspace)
print('Numero de puntos fijo:', equi_linspace) # N-1 puntos entre [0, 1] = {0, 1/n, 2/n, 3/n, ... , n/n = 1}
print('\n')
# Equiespaciados >> Elijo EL PASO entre los puntos en un dado intervalo [a, b)
equi_arange = np.arange(0, 1, 1./10) # como el range de Python pero puede venir con un paso en coma flotante
print('Paso entre puntos fijo:', equi_arange) # {0, paso, 2*paso, ... n*paso} siempre que n*paso<1
print('\n Identidad de 3×3')
identidad = np.identity(3)
print(identidad)
#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()
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 [9]:
x = np.linspace(0, 10, 1000) # ese array tiene 1000 elementos, andar printeando es poco práctico!
print(x.dtype) # Esto nos dice qué tipo de elemento es x (tener en cuenta que acá no van paréntesis)
ceros = np.zeros((100, 100)) # matriz de 100x100
print(ceros.shape) # np.shape nos dice cuántas filas y columnas tiene el array
# obs: esto es equivalente a:
# print(np.shape(ceros))
# Prueben qué pasa cuando le piden el shape a un array con una sola fila o columna como el x
La idea de este ejercicio es crear una matriz de 16x16 con los números del 1 al 256 siguiendo los pasos:
1) Crear un array que contenga los valores 1, 2, ..., 256
2) Modificar las dimensiones del array de 1×256 a 16×16
Ayuda: siempre que puedan utilicen las funciones de creación de arrays de numpy.
Ayuda bis: la función np.reshape hace lo que promete (reshapear) y puede tomar tuplas como argumento. Vean como usarla ya sea con su documentacion help(np.reshape) o con google.
A lo largo de nuestras carreras, y de la vida, nos encontramos con información que queremos graficar. 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)
In [10]:
import numpy as np
import matplotlib.pyplot as plt
# La siguiente linea va si estás usando un Notebook
%matplotlib inline
# 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() # Esta linea va si no estás usando un Notenook
Desglosemos el código que acabamos de escribir:
x y nuestra imagen y.plot de matplotlib.pyplot para graficar y en función de x.Aclaración 1: plt.show() vs %matplotlib:
En este ejemplo pusimos estas dos lineas, pero son redudantes, solo una de ellas es necesaria. Si están trabajando en un Notebook, o en algo que use un Notebook de fondo, este puede trabajar con gráficos de forma interactiva, por lo que necesitan de la linea %matplotlib inline. En el caso contrario, si están trabajando en un IDLE más usual, o ejecutando Python desde la terminal, necesitan usar la función plt.show() para mostrar los gráficos. No se vuelvan locos/as con esto, simplemente prueben cuál les funcina! (En Spyder puede funcionar cualquiera de las dos dependiendo de la configuración, pero por defecto se usa plt.show().)
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:
In [11]:
plt.plot(x, y, '.-');
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:
In [12]:
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!
Vemos que los puntos que se grafican son los $(x,y)$ = $\{(0,2);\,(2,4);\,(6,6);\,(8,14);\,(10,14)\}$
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:
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.
In [13]:
import numpy as np
import matplotlib.pyplot as plt
# La siguiente linea solo es necesaria si estás usando Jupyter Notebook
%matplotlib inline
# 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()
plt.show()
La mayoría de las funciones son bastante descriptivas. Veamos solo algunas aclaraciones:
plt.title es un argumento opcional y también se puede utilizar en el plt.xlabel o plt.ylabel por ejemplo.plt.grid() agrega una grid si no exite y la elimina si existe. Si se quiere forzar a que aparezca se puede utilizar plt.grid(True)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. (ver help(plt.legend))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í.
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('data_exp.csv', delimiter=',', skiprows=1, unpack=True)
Desglosemos esto por argumentos:
Por qué unpack=True
np.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. Utilizando unpack=True trasponemos 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 utilizando unpack=True y sin utilizarlo
También podríamos importarlo sin usar el unpack=True y luego usar data = np.transpose(data)
In [14]:
data = np.loadtxt('data_exp.csv', delimiter=',', skiprows=1)
print(f"Data sin unpack:\n {data}")
data = np.loadtxt('data_exp.csv', delimiter=',', skiprows=1, unpack=True)
print(f"\nData con unpack:\n {data}")
In [15]:
import numpy as np
import matplotlib.pyplot as plt
# Importamos los datos
data = np.loadtxt('data_exp.csv', delimiter=',', skiprows=1, unpack=True)
var_x, var_y = data # es lo mismo que poner en dos lineas x = data[0] e y = data[1]
# Graficamos los datos crudos
plt.plot(var_x, var_y, 'o-k', label='Datos')
plt.legend()
plt.show()
Los datos parecen tener un comportamiento exponencial. 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:
In [16]:
# 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(x, A, C): return A*np.exp(C*x)
# Utilizo curve_fit() para el ajuste
popt, pcov = curve_fit( f_ajuste, var_x, var_y )
# 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}')
# Graficamos los datos crudos
plt.plot(var_x, var_y, 'o-k', label='Datos')
# Graficamos el ajuste
A, C = popt
var_x_ajuste = var_x
var_y_ajuste = f_ajuste(var_x, A, C)
plt.plot(var_x_ajuste, var_y_ajuste, '-r', label='Ajuste')
# Legend y show
plt.legend()
plt.show()
El primer argumento es la función modelo con la que queremos ajustar.
La función modelo tiene que tener como primer argumento la variable que corresponde al eje $\hat x$. Luego se usan los otros argumentos para poner las constantes que queremos hallar (en nuestro caso A y C).
Los otros dos argumentos son nuestros datos
Estos son los datos crudos, los que usamos para visualizarlos normalmente.
curve_fit()popt
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 (A y C). Estos valores vienen ordenados del mismo modo que están en la función modelo por lo que popt[0] es el A óptimo y popt[1] es el C óptimo.
pcov
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ía
err_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, lo que implica que el dominio para graficar nuestro 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.
curve_fit() ven que como parámetro opcional también se le puede agregar el error de nuestros datos.plt.errorbar(), que funciona de forma muy similar a plt.plot() pero nos permite dar un array con los errores o un error constante para todos.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')
Perfecto, utilizando todo lo aprendido el código para gráfico completo y exportado sería:
In [17]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# Importamos los datos
var_x, var_y = np.loadtxt('data_exp.csv', delimiter=',', skiprows=1, unpack=True)
err_var_y = 1.1 * np.ones(len(var_y)) # ponemos que todos los var_y tiene 1.1 de error a modo de ej
# Graficamos los datos crudos con su error
plt.errorbar(var_x, var_y, yerr=err_var_y, fmt='o-k', label='Datos')
# Calculamos el ajuste
def f_ajuste(x, A, C): return A*np.exp(C*x)
popt, pcov = curve_fit(f_ajuste, var_x, var_y, sigma=err_var_y)
A, C = popt
# Declaramos nuestro nuevo dominio e imagen y graficamos el ajuste
new_var_x = np.linspace(min(var_x), max(var_x), 1000)
new_var_y = f_ajuste(new_var_x, A, C)
plt.plot(new_var_x, new_var_y, '-r', label='Ajuste exponencial')
# Títuo y labels
plt.title('Datos y Ajuste Exponencial', fontsize=16)
plt.xlabel('Variable 1', fontsize=14)
plt.ylabel('Variable 2', fontsize=14)
# Grid, legend, save y show
plt.grid(True)
plt.legend()
plt.savefig('ajuste_exponencial.png')
plt.show()
# Printiamos en pantalla los parámetros óptimos con sus errores
err_A, err_C = np.sqrt(np.diag(pcov))
print(f'A: {A} ± {err_A}')
print(f'C: {C} ± {err_C}')
La idea ahora es que hagan un script similar al de arriba para un set de datos lineal. Dejamos una pequeña guía.
Este ejercicio es el más largo y complicado hasta ahora. Hay muchas cosas! No se preocupen si no llegan, después en su casa pueden ir paso a paso leyendo la guía, probando y preguntando cosas por el canal de Discord.
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 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.
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.
In [18]:
import pandas as pd
df = pd.DataFrame(np.random.randn(8, 3), columns=['A', 'B', 'Columna'])
print(df) # Esto lo muestra mas "crudamente"
print()
df # Si solo escribimos esto la Notebook lo printea re cheto
# pero solo la notebook (Y quizas la terminal de spyder)
Out[18]:
Además podemos acceder a cada columna con su nombre! Genial. De yapa nos dice el tipo de dato que esa columna contiene.
In [19]:
print(df.A)
print()
print(df.Columna)
Ya que está podemos indexear de otra forma y no con números del 0-7. Abajo mostramos una forma para indexear con fechas consecutivas. Además podemos tomar cachos del DataFrame igual que como hacíamos con vectores antes.
In [20]:
dates = pd.date_range('20200404', periods=8)
df = pd.DataFrame({'Columna': np.random.randint(0,5,8),
'B': np.random.rand(8)},
index=dates)
print(df)
print()
print(df.loc[:, ['B']])
print()
print(df[0:4])
print()
print(df.sort_values(by='Columna'))
La librería también viene con una forma bastante poderosa para levantar archivos de un montón de tipos distintos (incluidas cosas tipo excel, json, html, etc...), vamos a usar read_csv para levantar el archivo que vimos arriba de la exponencial.
In [21]:
pd.read_csv("data_exp.csv")
#
# Escrito asi nomas dejamos que Pandas trate de leer el archivo lo mejor que pueda
# pero hay muchas opciones para levantar archivos mas "conflictivos", como para
# definir el separador entre datos, las lineas iniciales, tipo de datos, etc...
# esto puede verse desde la documentación de la librería o google.
#
Out[21]:
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 alrededor de algún valor medio usando randn 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.
In [22]:
import numpy as np
import matplotlib.pyplot as plt
mu = 100 # mu es mi valor medio
sigma = 15 # sigma la desviación
x = mu + sigma*np.random.randn(10000) # le sumo ruido gaussiano a mu
n, bins, patches = plt.hist(x, bins=50, edgecolor='black', facecolor='green', 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
In [23]:
# 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(x, bins=50, edgecolor='black', facecolor='green', alpha=0.75, density=True)
# Ver la escala vertical.
Y ya que estamos, para concientizar acerca de los peligros a la hora de la elección de bins, graficamos algunos histogramas superpuestos.
In [24]:
n, bins, patches = plt.hist(x, bins=100, density=True, edgecolor='black', facecolor='green', alpha=0.75)
n, bins, patches = plt.hist(x, bins=10 , density=True, edgecolor='black', facecolor='red' , alpha=0.5)
n, bins, patches = plt.hist(x, bins=2 , density=True, edgecolor='black', facecolor='yellow' , alpha=0.3)
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, y esa función utiliza el 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 def que tome como primer argumento $\vec{X}$, segundo argumento $t$ (por ser la variable independiente) y luego los demas parametros que necesite (en este caso sean $g$ y $l$) y opere para obtener $\dot{\vec{X}}$.
In [25]:
from scipy.integrate import odeint
import numpy as np
In [26]:
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 con la biblioteca ipywidgets o tantas otras posibles.
In [27]:
# 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
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:
In [28]:
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
Y para un sistema de ecuaciones del tipo $Ax = b$:
In [29]:
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
Por supuesto, también se puede hacer producto matriz con vector, y... oh si, se pueden calcular inversas.
In [30]:
print(np.dot(mat,x))
print(linalg.inv(mat))
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} $$
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?
Ahora no sólo nos vamos a emancipar del Origin, sino también del bendito Wolfram Alpha. 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.
In [31]:
from sympy 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 (para quienes lo hayan hecho).
Una forma conveniente para empezar es utilizando el siguiente comando:
In [32]:
init_session(use_latex='matplotlib') #el argumento nos permite renderear los resultados con el latex de matplotlib
Vimos aquí que una cantidad de variables se definieron como symbols de distintos tipos. Veamos ahora algunas de las posibilidades que tenemos:
In [33]:
expand( (x + y)**5 )
Out[33]:
In [34]:
factor( x**6 - 1 )
Out[34]:
Utilizando variables simbólicas de tipo integer, podemos hacer algunas sumatorias, por ejemplo, la famosa $$\sum_{k=0}^{m}k$$
In [35]:
Sum(k, (k,0,m) ).doit().factor() #el comando .doit() evalúa la suma
Out[35]:
O incluso algunas series, como $$ \sum_{n=1}^{\infty}\frac{1}{n^2} $$
In [36]:
Sum(1/n**2, (n, 1, oo)).doit() #el infinito se escribe como oo (dos o minúscula)
Out[36]:
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}$$
In [37]:
diff(x**n , x).simplify()
Out[37]:
In [38]:
diff(sin(tan(8*x)), x)
Out[38]:
In [39]:
diff(sin(tan(8*x)), x).simplify()
Out[39]:
In [40]:
diff(x*sin(y), x, y) #qué pasa si en lugar de x,y ponemos y,y?
Out[40]:
In [41]:
integrate(1/(1+x**2), x)
Out[41]:
In [42]:
integrate( sin(x)/x, (x,-oo,oo) )
Out[42]:
Calculen la siguiente integral: $$ \int_{-\infty}^{+\infty} e^{-x^2} \text{d}x $$ (opcional: verificar el resultado a mano!)
In [43]:
# Realicen el ejercicio 9
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
En este caso, la función a importar es quad, porque usa un método de cuadraturas
In [44]:
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$$
In [45]:
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$
In [46]:
a = -1
b = 1
integral, error = quad(integrando, a, b)
print(integral, error)
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$$
In [47]:
# 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
In [48]:
plt.figure()
plt.plot(x, integrando(x), 'b', label = 'integrando')
plt.plot(x, prim, 'r', label = 'primitiva')
plt.grid(True)
plt.legend()
Out[48]:
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
In [49]:
import scipy.special as sp
Empecemos graficando un ejemplo común en probabilidad, la función error:
In [50]:
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:
In [51]:
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')
Out[51]: