Primeros pasos con Jupyter Notebook


¿Qué es Jupyter Notebook?

Jupyter Notebook es una herramienta increíblemente poderosa para desarrollar y presentar proyectos de ciencia de datos de forma interactiva. Este post intentara guiarle de cómo usar Jupyter Notebooks para proyectos de ciencia de datos y cómo configurarlo en su máquina local.

jupyter-notebook las mejores herramientas gratuitas de ciencia de datos

Pero primero: ¿qué es un “cuaderno”?

Un cuaderno integra código y su salida en un solo documento que combina visualizaciones, texto narrativo, ecuaciones matemáticas y otros medios enriquecidos. En otras palabras: es un documento único en el que puede ejecutar código, mostrar el resultado y también agregar explicaciones, fórmulas, gráficos y hacer que su trabajo sea más transparente, comprensible, repetible y compartible. 

El uso de Notebooks es ahora una parte importante del flujo de trabajo de la ciencia de datos en empresas de todo el mundo. Si su objetivo es trabajar con datos, el uso de una computadora portátil acelerará su flujo de trabajo y facilitará la comunicación y el intercambio de resultados. 

Lo mejor de todo es que, como parte del proyecto de código abierto  Jupyter , los cuadernos de Jupyter son completamente gratuitos. Puede descargar el software  solo o como parte del kit de herramientas de ciencia de datos de Anaconda .

Aunque es posible utilizar muchos lenguajes de programación diferentes en Jupyter Notebooks, este artículo se centrará en Python, ya que es el caso de uso más común. (Entre los usuarios de R, R Studio  tiende a ser una opción más popular).Jupyter Notebooks también puede actuar como una plataforma flexible para familiarizarse con los pandas e incluso con Python, como se verá en este tutorial.

Algunas cosas que se veran en este post:

  • Cubrir los conceptos básicos de la instalación de Jupyter y la creación de su primer cuaderno
  • Profundizar y aprender toda la terminología importante
  • Explorear la facilidad con la que se pueden compartir y publicar en línea los blocs de notas.

Ejemplo de análisis de datos en un cuaderno Jupyter

Primero, recorreremos la configuración y un análisis de muestra para responder una pregunta de la vida real. Esto demostrará cómo el flujo de una computadora portátil hace que las tareas de ciencia de datos sean más intuitivas para nosotros mientras trabajamos y para otros una vez que llega el momento de compartir nuestro trabajo.

Entonces, digamos que es analista de datos y se le ha encomendado la tarea de averiguar cómo cambiaron históricamente las ganancias de las empresas más grandes de los EE. UU. Encontrará un conjunto de datos de compañías Fortune 500 que abarcan más de 50 años desde la primera publicación de la lista en 1955, reunidos a partir del archivo público de  Fortune . Seguimos adelante y creamos un CSV de los datos que puede usar  aquí .

Como demostraremos, los cuadernos de Jupyter son perfectamente adecuados para esta investigación. Primero, sigamos e instalemos Jupyter.

Instalación

La forma más fácil para que un principiante comience con Jupyter Notebooks es instalando  Anaconda .

Anaconda es la distribución de Python más utilizada para la ciencia de datos y viene precargada con todas las bibliotecas y herramientas más populares. 

Algunas de las bibliotecas de Python más grandes incluidas en Anaconda incluyen NumPy ,  pandas y  Matplotlib , aunque la  lista completa de más de 1000 es exhaustiva.

Por lo tanto, Anaconda nos permite comenzar con un taller de ciencia de datos completamente equipado sin la molestia de administrar innumerables instalaciones o preocuparse por las dependencias y los problemas de instalación específicos del sistema operativo (léase: específicos de Windows).

Para obtener Anaconda, simplemente:

  1. Descargue la última versión de Anaconda para Python 3.8.
  2. Instale Anaconda siguiendo las instrucciones en la página de descarga y / o en el ejecutable.

Si es un usuario más avanzado con Python ya instalado y prefiere administrar sus paquetes manualmente, puede usar pip :

pip3 install jupyter

Creación de su primer cuaderno

En esta sección, aprenderemos a ejecutar y guardar blocs de notas, a familiarizarnos con su estructura y a comprender la interfaz. Nos familiarizaremos con alguna terminología básica que lo guiará hacia una comprensión práctica de cómo usar Jupyter Notebooks por su cuenta y nos preparará para la siguiente sección, que recorre un ejemplo de análisis de datos y da vida a todo lo que aprendemos aquí.

Ejecutando Jupyter

En Windows, puede ejecutar Jupyter a través del acceso directo que Anaconda agrega a su menú de inicio, que abrirá una nueva pestaña en su navegador web predeterminado que debería parecerse a la siguiente captura de pantalla.

Panel de control de Jupyter

Este no es un cuaderno por el momento, ¡pero que no cunda el pánico! No hay mucho que hacer. Este es el panel de control del portátil, diseñado específicamente para administrar sus portátiles Jupyter. Piense en ello como la plataforma de lanzamiento para explorar, editar y crear sus cuadernos.

Tenga en cuenta que el panel le dará acceso solo a los archivos y subcarpetas que se encuentran dentro del directorio de inicio de Jupyter (es decir, donde está instalado Jupyter o Anaconda). Sin embargo, el directorio de inicio se puede cambiar .

También es posible iniciar el tablero en cualquier sistema a través del símbolo del sistema (o terminal en sistemas Unix) ingresando el comando jupyter notebook; en este caso, el directorio de trabajo actual será el directorio de inicio.

Con Jupyter Notebook abierto en su navegador, es posible que haya notado que la URL del panel de control es algo así como https://localhost:8888/tree. Localhost no es un sitio web, pero indica que el contenido se sirve desde su   máquina local : su propia computadora.

Los cuadernos y el tablero de Jupyter son aplicaciones web, y Jupyter inicia un servidor Python local para servir estas aplicaciones en su navegador web, lo que lo hace esencialmente independiente de la plataforma y abre la puerta para compartir más fácilmente en la web.

(Si aún no comprende esto, no se preocupe; el punto importante es que, aunque Jupyter Notebooks se abre en su navegador, está alojado y se ejecuta en su máquina local. Sus blocs de notas no están realmente en la web hasta que decide compartirlos.)

La interfaz del tablero se explica por sí misma, aunque volveremos a ella brevemente más adelante. entonces que estamos esperando ‘ Busque la carpeta en la que le gustaría crear su primer cuaderno, haga clic en el botón desplegable “Nuevo” en la parte superior derecha y seleccione “Python 3”:

Nuevo menú de cuaderno

¡Oye presto, aquí estamos! Su primer cuaderno de Jupyter se abrirá en una nueva pestaña; cada cuaderno usa su propia pestaña porque puede abrir varios cuadernos simultáneamente.

Si regresa al tablero, verá el nuevo archivo Untitled.ipynb y debería ver un texto verde que le indica que su computadora portátil se está ejecutando.

¿Qué es un archivo ipynb?

La respuesta corta: cada  archivo .ipynb es un cuaderno, por lo que cada vez que cree un nuevo cuaderno, se creará un nuevo   archivo .ipynb.  

La respuesta más larga: cada .ipynb archivo es un archivo de texto que describe el contenido de su cuaderno en un formato llamado  JSON . Cada celda y su contenido, incluidos los archivos adjuntos de imágenes que se han convertido en cadenas de texto, se enumeran allí junto con algunos  metadatos .

Puede editarlo usted mismo, ¡si sabe lo que está haciendo! – seleccionando “Editar> Editar metadatos del cuaderno” en la barra de menú del cuaderno. También puede ver el contenido de los archivos de su cuaderno seleccionando “Editar” en los controles del panel.

Sin embargo, la palabra clave puede. En la mayoría de los casos, no hay ninguna razón por la que deba editar manualmente los metadatos de su cuaderno.

La interfaz del portátil

Ahora que tiene un cuaderno abierto frente a usted, es de esperar que su interfaz no se vea del todo extraña. Después de todo, Jupyter es esencialmente un procesador de texto avanzado.

¿Por qué no echar un vistazo? Consulte los menús para familiarizarse con ellos, especialmente tómese unos minutos para desplazarse hacia abajo en la lista de comandos en la paleta de comandos, que es el botón pequeño con el ícono del teclado (o Ctrl + Shift + P).

Nuevo cuaderno de Jupyter

Hay dos términos bastante prominentes que debe notar, que probablemente sean nuevos para usted: las  células  y los  núcleos  son clave tanto para comprender Jupyter como para lo que lo convierte en algo más que un procesador de texto. Afortunadamente, estos conceptos no son difíciles de entender.

  • Un kernel es un “motor computacional” que ejecuta el código contenido en un documento de cuaderno.
  • Una celda es un contenedor para el texto que se mostrará en el cuaderno o el código que ejecutará el kernel del cuaderno.

Celdas

Regresaremos a los núcleos un poco más tarde, pero primero hablemos de las celdas. Las células forman el cuerpo de un cuaderno. En la captura de pantalla de un nuevo cuaderno en la sección anterior, ese cuadro con el contorno verde es una celda vacía. Hay dos tipos de células principales que cubriremos:

  • Una  celda de código contiene código que se ejecutará en el kernel. Cuando se ejecuta el código, el portátil muestra el resultado debajo de la celda de código que lo generó.
  • Una  celda de Markdown contiene texto formateado con Markdown y muestra su salida en el lugar cuando se ejecuta la celda de Markdown.

La primera celda de un cuaderno nuevo es siempre una celda de código.

Probémoslo con un ejemplo clásico de Hello World: escriba print('Hello World!') en la celda y haga clic en el botón Ejecutar  Botón de ejecución del cuaderno en la barra de herramientas de arriba o presione  Ctrl + Enter.

El resultado debería verse así:

print('Hello World!')
Hello World!

Cuando ejecutamos la celda, su salida se muestra a continuación y la etiqueta a su izquierda habrá cambiado de In [ ] a  In [1].

La salida de una celda de código también forma parte del documento, por lo que puede verla en este artículo. Siempre puede notar la diferencia entre el código y las celdas de Markdown porque las celdas de código tienen esa etiqueta a la izquierda y las celdas de Markdown no.

La parte “En” de la etiqueta es simplemente la abreviatura de “Entrada”, mientras que el número de etiqueta indica cuándo se ejecutó la celda en el kernel; en este caso, la celda se ejecutó primero.

Ejecute la celda nuevamente y la etiqueta cambiará a In [2] porque ahora la celda fue la segunda en ejecutarse en el kernel. Será más claro por qué esto es tan útil más adelante cuando echemos un vistazo más de cerca a los núcleos.

En la barra de menú, haga clic en  Insertar  y seleccione  Insertar celda debajo  para crear una nueva celda de código debajo de la primera y pruebe el siguiente código para ver qué sucede. ¿Notas algo diferente?

import time
time.sleep(3)

Esta celda no produce ningún resultado, pero tarda tres segundos en ejecutarse. Observe cómo Jupyter significa cuando la celda se está ejecutando actualmente cambiando su etiqueta a In [*]

En general, la salida de una celda proviene de cualquier dato de texto impreso específicamente durante la ejecución de la celda, así como del valor de la última línea de la celda, ya sea una variable solitaria, una llamada de función u otra cosa. Por ejemplo:

def say_hello(recipient):
    return 'Hello, {}!'.format(recipient)

say_hello('Tim')
'Hello, Tim!'

Te encontrarás usando esto casi constantemente en tus propios proyectos, y veremos más de eso más adelante.

Atajos de teclado

Una última cosa que puede haber observado al ejecutar sus celdas es que su borde se vuelve azul, mientras que era verde mientras estaba editando. En un Jupyter Notebook, siempre hay una celda “activa” resaltada con un borde cuyo color denota su modo actual:

  • Contorno verde : la celda está en “modo de edición”
  • Contorno azul : la celda está en “modo de comando”

Entonces, ¿qué podemos hacer con una celda cuando está en modo comando? Hasta ahora, hemos visto cómo ejecutar una celda  Ctrl + Enter, pero hay muchos otros comandos que podemos usar. La mejor manera de usarlos es con atajos de teclado.

Los atajos de teclado son un aspecto muy popular del entorno de Jupyter porque facilitan un flujo de trabajo rápido basado en celdas. Muchas de estas son acciones que puede realizar en la celda activa cuando está en modo de comando.

A continuación, encontrará una lista de algunos de los atajos de teclado de Jupyter. No es necesario que los memorice todos de inmediato, pero esta lista debería darle una buena idea de lo que es posible.

  • Alternar entre los modos de edición y comando con  Esc y  Enter, respectivamente.
  • Una vez en el modo de comando:
    • Desplácese hacia arriba y hacia abajo en sus celdas con las   teclas Up y  Down.
    • Presione  A o  B para insertar una nueva celda encima o debajo de la celda activa.
    • M transformará la celda activa en una celda de Markdown.
    • Y establecerá la celda activa en una celda de código.
    • D + D ( D dos veces) eliminará la celda activa.
    • Z deshará la eliminación de la celda.
    • Mantenga  Shift presionado y presione  Up o  Downpara seleccionar varias celdas a la vez. Con varias celdas seleccionadas, Shift + M fusionará su selección.
  • Ctrl + Shift + -, en el modo de edición, dividirá la celda activa en el cursor.
  • También puede hacer clic Shift + Click en y  en el margen a la izquierda de sus celdas para seleccionarlas.

Continúe y pruébelos en su propio cuaderno. Una vez que esté listo, cree una nueva celda de Markdown y aprenderemos a formatear el texto en nuestros cuadernos.

Reducción

Markdown es un lenguaje de marcado ligero y fácil de aprender para formatear texto sin formato. Su sintaxis tiene una correspondencia uno a uno con las etiquetas HTML, por lo que algunos conocimientos previos aquí serían útiles, pero definitivamente no son un requisito previo.

Recuerda que este artículo fue escrito en un cuaderno de Jupyter, por lo que todo el texto narrativo y las imágenes que has visto hasta ahora se lograron escribiendo en Markdown. Cubramos los conceptos básicos con un ejemplo rápido:

# This is a level 1 heading

## This is a level 2 heading

This is some plain text that forms a paragraph. Add emphasis via **bold** and __bold__, or *italic* and _italic_. 

Paragraphs must be separated by an empty line. 

* Sometimes we want to include lists. 
* Which can be bulleted using asterisks. 

1. Lists can also be numbered. 
2. If we want an ordered list.

[It is possible to include hyperlinks](https://www.example.com)

Inline code uses single backticks: `foo()`, and code blocks use triple backticks: 
```
bar()
``` 
Or can be indented by 4 spaces: 

    foo()
    
And finally, adding images is easy: ![Alt text](https://www.example.com/image.jpg)

Así es como se vería Markdown una vez que ejecute la celda para renderizarlo:

(Tenga en cuenta que el texto alternativo de la imagen se muestra aquí porque en realidad no usamos una URL de imagen válida en nuestro ejemplo)

Al adjuntar imágenes, tiene tres opciones:

  • Utilice una URL para una imagen en la web.
  • Use una URL local para una imagen que mantendrá junto a su cuaderno, como en el mismo repositorio de git.
  • Agregue un archivo adjunto a través de “Editar> Insertar imagen”; esto convertirá la imagen en una cadena y la almacenará dentro de su .ipynbarchivo de cuaderno  . ¡Tenga en cuenta que esto hará que su .ipynb archivo sea mucho más grande!

Hay mucho más en Markdown, especialmente en torno a los hipervínculos, y también es posible incluir simplemente HTML simple. Una vez que se encuentre superando los límites de los conceptos básicos anteriores, puede consultar la guía oficial  del creador de Markdown, John Gruber, en su sitio web.

Granos

Detrás de cada portátil hay un kernel. Cuando ejecuta una celda de código, ese código se ejecuta dentro del kernel. Cualquier salida se devuelve a la celda para que se muestre. El estado del kernel persiste a lo largo del tiempo y entre celdas; pertenece al documento como un todo y no a celdas individuales.

Por ejemplo, si importa bibliotecas o declara variables en una celda, estarán disponibles en otra. Probemos esto para sentirlo. Primero, importaremos un paquete de Python y definiremos una función:

import numpy as np
def square(x):
    return x * x

Una vez que hayamos ejecutado la celda de arriba, podemos hacer referencia a np y  squareen cualquier otra celda. 

x = np.random.randint(1, 10)
y = square(x)
print('%d squared is %d' % (x, y))
1 squared is 1

Esto funcionará independientemente del orden de las celdas de su cuaderno. Siempre que se haya ejecutado una celda, las variables que declaró o las bibliotecas que importó estarán disponibles en otras celdas.

Puede probarlo usted mismo, imprimamos nuestras variables nuevamente.

print('Is %d squared %d?' % (x, y))
Is 1 squared 1?

¡No hay sorpresas aquí! Pero, ¿qué pasa si cambiamos el valor de  y?

y = 10
print('Is %d squared is %d?' % (x, y))

Si ejecutamos la celda de arriba, ¿qué crees que pasaría?

Obtendremos una salida como:  Is 4 squared 10?. Esto se debe a que una vez que hemos ejecutado la  y = 10celda de código, yya no es igual al cuadrado de x en el kernel. 

La mayoría de las veces, cuando crea un cuaderno, el flujo será de arriba hacia abajo. Pero es común volver atrás para realizar cambios. Cuando necesitamos realizar cambios en una celda anterior, el orden de ejecución que podemos ver a la izquierda de cada celda, por ejemplo In [6], puede ayudarnos a diagnosticar problemas al ver en qué orden se han ejecutado las celdas. 

Y si alguna vez deseamos restablecer las cosas, hay varias opciones increíblemente útiles en el menú Kernel:

  • Reiniciar: reinicia el kernel, borrando así todas las variables, etc.que fueron definidas.
  • Reiniciar y borrar la salida: igual que el anterior, pero también borrará la salida que se muestra debajo de las celdas de código.
  • Reiniciar y ejecutar todo: igual que el anterior, pero también ejecutará todas sus celdas en orden de la primera a la última.

Si su kernel se atasca alguna vez en un cálculo y desea detenerlo, puede elegir la opción Interrumpir.

Elegir un kernel

Es posible que haya notado que Jupyter le brinda la opción de cambiar el kernel y, de hecho, hay muchas opciones diferentes para elegir. Cuando creó un nuevo cuaderno desde el panel de control seleccionando una versión de Python, en realidad estaba eligiendo qué kernel usar.

Hay kernels para diferentes versiones de Python, y también para más de 100 lenguajes,  incluidos Java, C e incluso Fortran. Los científicos de datos pueden estar particularmente interesados ​​en los núcleos para  R  y  Julia , así como en  imatlab  y  Calysto MATLAB Kernel  para Matlab.

El kernel de SoS  proporciona compatibilidad con varios idiomas en un solo portátil.

Cada kernel tiene sus propias instrucciones de instalación, pero probablemente requiera que ejecute algunos comandos en su computadora.

¿Aprendiendo Python para la ciencia de datos?

¡Prueba una de nuestras más de 60 misiones gratuitas hoy!Empiece >>

Análisis de ejemplo

Ahora que hemos visto  qué es  un Jupyter Notebook, es hora de ver  cómo se utilizan en la práctica, lo que debería darnos una comprensión más clara de por  qué son tan populares.

Finalmente es hora de comenzar con ese conjunto de datos de Fortune 500 mencionado anteriormente. Recuerde, nuestro objetivo es averiguar cómo han cambiado históricamente las ganancias de las empresas más grandes de EE . UU .

Vale la pena señalar que todos desarrollarán sus propias preferencias y estilo, pero los principios generales aún se aplican. Si lo desea, puede seguir esta sección en su propio cuaderno o utilizarla como guía para crear su propio enfoque.

Nombrar sus cuadernos

Antes de comenzar a escribir su proyecto, probablemente querrá darle un nombre significativo. nombre de archivo Untitleden la parte superior izquierda de la pantalla para ingresar un nuevo nombre de archivo, y presione el icono Guardar (que parece un disquete) debajo de él para guardar.

Tenga en cuenta que cerrar la pestaña del cuaderno en su navegador no “cerrará” su cuaderno de la forma en que lo haría cerrar un documento en una aplicación tradicional. El kernel del portátil seguirá ejecutándose en segundo plano y debe cerrarse antes de que se “cierre” realmente, aunque esto es muy útil si cierra accidentalmente la pestaña o el navegador.

Si el kernel se apaga, puede cerrar la pestaña sin preocuparse de si todavía se está ejecutando o no. 

La forma más sencilla de hacerlo es seleccionar “Archivo> Cerrar y detener” en el menú de la libreta. Sin embargo, también puede apagar el kernel yendo a “Kernel> Shutdown” desde la aplicación de la computadora portátil o seleccionando la computadora portátil en el panel de control y haciendo clic en “Apagar” (vea la imagen a continuación).

Un cuaderno para correr

Configuración

Es común comenzar con una celda de código específicamente para las importaciones y la configuración, de modo que si elige agregar o cambiar algo, simplemente puede editar y volver a ejecutar la celda sin causar efectos secundarios.

%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns sns.set(style="darkgrid")

Importaremos pandas  para trabajar con nuestros datos,  Matplotlib  para trazar gráficos y  Seaborn  para hacer que nuestros gráficos sean más bonitos. También es común importar  NumPy, pero en este caso, pandas lo importa por nosotros.

Esa primera línea no es un comando de Python, pero usa algo llamado magia de línea para indicarle a Jupyter que capture los diagramas de Matplotlib y los represente en la salida de la celda. Hablaremos un poco más sobre la magia de líneas más adelante, y también se tratan en nuestro tutorial avanzado de Jupyter Notebooks .

Por ahora, sigamos adelante y carguemos nuestros datos.

df = pd.read_csv('fortune500.csv')

Es sensato hacer esto también en una sola celda, en caso de que necesitemos recargarlo en algún momento.

Guardar y comprobar

Ahora que hemos comenzado, lo mejor es ahorrar con regularidad. Al presionar  Ctrl + Sguardará nuestro cuaderno llamando al comando “Guardar y punto de control”, pero ¿qué es este punto de control?

Cada vez que creamos un nuevo cuaderno, se crea un archivo de punto de control junto con el archivo del cuaderno. Se encuentra dentro de un subdirectorio oculto de su ubicación de guardado llamado .ipynb_checkpoints y también es un  .ipynb archivo.

De forma predeterminada, Jupyter guardará automáticamente su cuaderno cada 120 segundos en este archivo de punto de control sin alterar el archivo de su cuaderno principal. Cuando “Guardar y punto de control”, se actualizan tanto el cuaderno como los archivos del punto de control. Por lo tanto, el punto de control le permite recuperar su trabajo no guardado en caso de un problema inesperado.

Puede volver al punto de control desde el menú a través de “Archivo> Volver al punto de control”.

Investigando nuestro conjunto de datos

¡Ahora estamos realmente rodando! Nuestro cuaderno se guarda de forma segura y hemos cargado nuestro conjunto de datos  df en la estructura de datos de pandas más utilizada, que se llama  DataFrame ay básicamente se parece a una tabla. ¿Qué aspecto tiene el nuestro?

df.head()
AñoRangoEmpresaIngresos (en millones)Beneficio (en millones)
019551Motores generales9823.5806
119552Exxon Mobil5661.4584,8
219553Acero de EE. UU.3250,4195,4
319554Energia General2959,1212,6
419555Esmark2510,819,1
df.tail()
AñoRangoEmpresaIngresos (en millones)Beneficio (en millones)
254952005496Wm. Wrigley Jr.3648,6493
254962005497Energía Peabody3631,6175,4
254972005498Wendy’s internacional3630,457,8
254982005499Kindred Healthcare3616,670,6
254992005500Cincinnati Financial3614.0584

Luciendo bien. Tenemos las columnas que necesitamos y cada fila corresponde a una sola empresa en un solo año.

Cambiemos el nombre de esas columnas para poder consultarlas más adelante.

df.columns = ['year', 'rank', 'company', 'revenue', 'profit']

A continuación, necesitamos explorar nuestro conjunto de datos. Esta completo? ¿Lo leyeron los pandas como se esperaba? ¿Falta algún valor?

len(df)
25500

De acuerdo, eso se ve bien, son 500 filas por cada año desde 1955 hasta 2005, inclusive.

Comprobemos si nuestro conjunto de datos se ha importado como era de esperar. Una simple verificación es ver si los tipos de datos (o dtypes) se han interpretado correctamente.

df.dtypes
year int64 rank int64 company object revenue float64 profit object dtype: object

UH oh. Parece que hay algún problema con la columna de ganancias; esperaríamos que fuera  float64 similar a la columna de ingresos. Esto indica que probablemente contiene algunos valores que no son enteros, así que echemos un vistazo.

non_numberic_profits = df.profit.str.contains('[^0-9.-]')
df.loc[non_numberic_profits].head()
añorangoempresaingresoslucro
2281955229Norton135,0N / A
2901955291Elaboración de cerveza Schlitz100,0N / A
2941955295Aceite vegetal del Pacífico97,9N / A
2961955297Cervecerías Liebmann96,0N / A
3521955353Minneapolis-Moline77,4N / A

¡Tal como sospechábamos! Algunos de los valores son cadenas, que se han utilizado para indicar datos faltantes. ¿Hay otros valores que se hayan infiltrado?

set(df.profit[non_numberic_profits])
{'N.A.'}

Eso hace que sea fácil de interpretar, pero ¿qué debemos hacer? Bueno, eso depende de cuántos valores falten.

len(df.profit[non_numberic_profits])
369

Es una pequeña fracción de nuestro conjunto de datos, aunque no del todo intrascendente, ya que todavía está en torno al 1,5%.

Si las filas que contienen N.A. están distribuidas de manera aproximada y uniforme a lo largo de los años, la solución más sencilla sería eliminarlas. Así que echemos un vistazo rápido a la distribución.

bin_sizes, _, _ = plt.hist(df.year[non_numberic_profits], bins=range(1955, 2006))
Falta distribución de valor

De un vistazo, podemos ver que los valores más inválidos en un solo año son menos de 25, y como hay 500 puntos de datos por año, eliminar estos valores representaría menos del 4% de los datos de los peores años. De hecho, aparte de un aumento en torno a los 90, la mayoría de los años tienen menos de la mitad de los valores faltantes del pico.

Para nuestros propósitos, digamos que esto es aceptable y continúe y elimine estas filas.

df = df.loc[~non_numberic_profits]
df.profit = df.profit.apply(pd.to_numeric)

Deberíamos comprobar que funcionó.

len(df)
25131
df.dtypes
year int64 rank int64 company object revenue float64 profit float64 dtype: object

¡Estupendo! Hemos terminado la configuración de nuestro conjunto de datos.

Si tuviéramos que presentar su cuaderno como un informe, podríamos deshacernos de las celdas de investigación que creamos, que se incluyen aquí como una demostración del flujo de trabajo con los cuadernos, y fusionar las celdas relevantes (consulte la sección Funcionalidad avanzada a continuación para más sobre esto) para crear una celda de configuración de un solo conjunto de datos.

Esto significaría que si alguna vez estropeamos nuestro conjunto de datos en otro lugar, podemos simplemente volver a ejecutar la celda de configuración para restaurarlo.

Trazando con matplotlib

A continuación, podemos llegar a abordar la cuestión en cuestión trazando el beneficio promedio por año. También podríamos trazar los ingresos, así que primero podemos definir algunas variables y un método para reducir nuestro código.

group_by_year = df.loc[:, ['year', 'revenue', 'profit']].groupby('year')
avgs = group_by_year.mean()
x = avgs.index
y1 = avgs.profit
def plot(x, y, ax, title, y_label):
    ax.set_title(title)
    ax.set_ylabel(y_label)
    ax.plot(x, y)
    ax.margins(x=0, y=0)

¡Ahora tramemos!

fig, ax = plt.subplots()
plot(x, y1, ax, 'Increase in mean Fortune 500 company profits from 1955 to 2005', 'Profit (millions)')
Aumento de las ganancias medias de las empresas incluidas en la lista Fortune 500 de 1955 a 2005

Vaya, eso parece exponencial, pero tiene algunas caídas enormes. Deben corresponder a la recesión de principios de la  década de 1990  y a la  burbuja de las puntocom . Es bastante interesante ver eso en los datos. Pero, ¿cómo es que las ganancias se recuperaron a niveles aún más altos después de cada recesión?

Quizás los ingresos puedan decirnos más.

y2 = avgs.revenue
fig, ax = plt.subplots()
plot(x, y2, ax, 'Increase in mean Fortune 500 company revenues from 1955 to 2005', 'Revenue (millions)')
Aumento de los ingresos medios de las empresas incluidas en la lista Fortune 500 de 1955 a 2005

Eso agrega otro lado a la historia. Los ingresos no se vieron tan afectados; ese es un gran trabajo contable de los departamentos de finanzas.

Con un poco de ayuda  de Stack Overflow , podemos superponer estos gráficos con +/- sus desviaciones estándar.

def plot_with_std(x, y, stds, ax, title, y_label):
    ax.fill_between(x, y - stds, y + stds, alpha=0.2)
    plot(x, y, ax, title, y_label)
fig, (ax1, ax2) = plt.subplots(ncols=2)
title = 'Increase in mean and std Fortune 500 company %s from 1955 to 2005'
stds1 = group_by_year.std().profit.values
stds2 = group_by_year.std().revenue.values
plot_with_std(x, y1.values, stds1, ax1, title % 'profits', 'Profit (millions)')
plot_with_std(x, y2.values, stds2, ax2, title % 'revenues', 'Revenue (millions)')
fig.set_size_inches(14, 4)
fig.tight_layout()
jupyter-notebook-tutorial_48_0

Eso es asombroso, ¡las desviaciones estándar son enormes! Algunas empresas de Fortune 500 ganan miles de millones mientras que otras pierden miles de millones, y el riesgo ha aumentado junto con el aumento de las ganancias a lo largo de los años.

Quizás algunas empresas se desempeñen mejor que otras; ¿Son las ganancias del 10% superior más o menos volátiles que las del 10% inferior?

Hay muchas preguntas que podríamos analizar a continuación, y es fácil ver cómo el flujo de trabajo en un cuaderno puede coincidir con el propio proceso de pensamiento. Para los propósitos de este tutorial, detendremos nuestro análisis aquí, ¡pero siéntase libre de continuar investigando los datos por su cuenta!

Este flujo nos ayudó a investigar fácilmente nuestro conjunto de datos en un solo lugar sin cambiar de contexto entre aplicaciones, y nuestro trabajo se puede compartir y reproducir de inmediato. Si quisiéramos crear un informe más conciso para una audiencia en particular, podríamos refactorizar rápidamente nuestro trabajo fusionando celdas y eliminando el código intermediario.

Compartir sus cuadernos

Cuando la gente habla de compartir sus cuadernos, generalmente hay dos paradigmas que pueden estar considerando.

La mayoría de las veces, las personas comparten el resultado final de su trabajo, al igual que este artículo, lo que significa compartir versiones no interactivas y pre-renderizadas de sus cuadernos. Sin embargo, también es posible colaborar en portátiles con la ayuda de sistemas de control de versiones como Git o plataformas online como Google Colab .

Antes de compartir

Un cuaderno compartido aparecerá exactamente en el estado en el que se encontraba cuando lo exporta o lo guarda, incluida la salida de las celdas de código. Por lo tanto, para asegurarse de que su computadora portátil esté lista para compartir, por así decirlo, hay algunos pasos que debe seguir antes de compartir:

  1. Haga clic en “Celda> Todos los resultados> Borrar”.
  2. Haga clic en “Kernel> Reiniciar y ejecutar todo”.
  3. Espere a que las celdas de su código terminen de ejecutarse y verifique que se haya ejecutado como se esperaba

Esto asegurará que sus cuadernos no contengan salida intermedia, tengan un estado obsoleto y se ejecuten en orden en el momento de compartir.

Exportación de sus cuadernos

Jupyter tiene soporte integrado para exportar a HTML y PDF, así como a varios otros formatos, que puede encontrar en el menú en “Archivo> Descargar como”.

Si desea compartir sus cuadernos con un pequeño grupo privado, esta funcionalidad puede ser todo lo que necesite. De hecho, como a muchos investigadores de instituciones académicas se les proporciona un espacio web público o interno, y debido a que puede exportar un cuaderno a un archivo HTML, Jupyter Notebooks puede ser una forma especialmente conveniente para que los investigadores compartan sus resultados con sus pares.

Pero si compartir archivos exportados no es suficiente para usted, también existen algunos métodos inmensamente populares para compartir  .ipynb archivos más directamente en la web.

GitHub

Con la  cantidad de cuadernos públicos en GitHub que  excedieron los 1.8 millones a principios de 2018, seguramente es la plataforma independiente más popular para compartir proyectos de Jupyter con el mundo. GitHub tiene soporte integrado para renderizar  .ipynb archivos directamente tanto en repositorios como en esencias en su sitio web. Si aún no lo sabe,  GitHub  es una plataforma de alojamiento de código para el control de versiones y la colaboración para repositorios creados con  Git. . Necesitará una cuenta para utilizar sus servicios, pero las cuentas estándar son gratuitas.

Una vez que tenga una cuenta de GitHub, la forma más fácil de compartir un cuaderno en GitHub no requiere Git en absoluto. Desde 2008, GitHub ha proporcionado su servicio Gist para alojar y compartir fragmentos de código, cada uno de los cuales tiene su propio repositorio. Para compartir un cuaderno usando Gists:

  1. Inicie sesión y navegue hasta gist.github.com .
  2. Abra su  .ipynb archivo en un editor de texto, seleccione todo y copie el JSON dentro.
  3. Pegue el JSON del cuaderno en la esencia.
  4. Dale a tu Gist un nombre de archivo, recordando agregarlo  .iypnb o esto no funcionará.
  5. Haga clic en “Crear esencia secreta” o “Crear esencia pública”.

Esto debería tener un aspecto similar al siguiente:

Creando una esencia

Si creó un Gist público, ahora podrá compartir su URL con cualquier persona, y otros podrán  bifurcar y clonar  su trabajo.

Crear tu propio repositorio de Git y compartirlo en GitHub está más allá del alcance de este tutorial, pero  GitHub proporciona muchas guías  para que comiences por tu cuenta.

Un consejo adicional para aquellos que usan git es  agregar una excepción  a su  .gitignore para los .ipynb_checkpoints directorios ocultos que  crea Jupyter, para no enviar archivos de puntos de control innecesariamente a su repositorio.

Nbviewer

Después de haber crecido hasta renderizar  cientos de miles  de cuadernos cada semana en 2015, NBViewer es el renderizador de cuadernos más popular en la web. Si ya tiene un lugar para alojar sus cuadernos de Jupyter en línea, ya sea en GitHub o en otro lugar, NBViewer renderizará su cuaderno y proporcionará una URL para compartir junto con él. Se proporciona como un servicio gratuito como parte del Proyecto Jupyter y está disponible en  nbviewer.jupyter.org .

Desarrollado inicialmente antes de la integración de Jupyter Notebook de GitHub, NBViewer permite que cualquier persona ingrese una URL, ID de Gist o nombre de usuario / repositorio / archivo de GitHub y mostrará el cuaderno como una página web. La ID de un Gist es el número único al final de su URL; por ejemplo, la cadena de caracteres después de la última barra invertida  https://gist.github.com/username/50896401c23e0bf417e89cd57e89e1de. Si ingresa un nombre de usuario o nombre de usuario / repositorio de GitHub, verá un explorador de archivos mínimo que le permite explorar los repositorios de un usuario y su contenido.

La URL que se muestra en NBViewer cuando se muestra un cuaderno es una constante basada en la URL del cuaderno que se está procesando, por lo que puede compartirla con cualquier persona y funcionará siempre que los archivos originales permanezcan en línea; NBViewer no almacena en caché los archivos durante mucho tiempo. largo.

Si no le gusta Nbviewer, existen otras opciones similares: aquí hay un hilo con algunos para considerar de nuestra comunidad.

Extras: Extensiones de Jupyter Notebook

Ya hemos cubierto todo lo que necesita para empezar a trabajar en Jupyter Notebooks.

¿Qué son las extensiones?

Las extensiones son precisamente lo que parecen: características adicionales que amplían la funcionalidad de Jupyter Notebooks. Si bien un Jupyter Notebook básico puede hacer mucho, las extensiones ofrecen algunas características adicionales que pueden ayudar con flujos de trabajo específicos o que simplemente mejoran la experiencia del usuario.

Por ejemplo, una extensión llamada “Tabla de contenido” genera una tabla de contenido para su cuaderno, para que los cuadernos grandes sean más fáciles de visualizar y navegar. 

Otro, llamado Inspector de variables, le mostrará el valor, el tipo, el tamaño y la forma de cada variable en su cuaderno para una fácil referencia rápida y depuración. 

Otro, llamado ExecuteTime, le permite saber cuándo y durante cuánto tiempo se ejecutó cada celda; esto puede ser particularmente conveniente si está tratando de acelerar un fragmento de su código.

Estos son solo la punta del iceberg; hay muchas extensiones disponibles.

¿Dónde puede obtener extensiones?

Para obtener las extensiones, debe instalar Nbextensions. Puede hacer esto usando pip y la línea de comando. Si tiene Anaconda, puede ser mejor hacerlo a través de Anaconda Prompt en lugar de la línea de comandos normal.

Cerrar Jupyter cuadernos, abierta Anaconda Prompt y ejecute el siguiente comando: pip install jupyter_contrib_nbextensions && jupyter contrib nbextension install.

Una vez que haya hecho eso, inicie un cuaderno y debería ver una pestaña Nbextensions. Al hacer clic en esta pestaña, se le mostrará una lista de extensiones disponibles. Simplemente marque las casillas de las extensiones que desea habilitar y ¡listo para las carreras!

Instalación de extensiones

Una vez que se ha instalado Nbextensions, no es necesario realizar una instalación adicional de cada extensión. Sin embargo, si ya instaló Nbextensons pero no ve la pestaña, no está solo. Este hilo en Github detalla algunos problemas y soluciones comunes.

Extras: Line Magics en Jupyter

Anteriormente mencionamos los comandos mágicos cuando solíamos %matplotlib inlinehacer que los gráficos de Matplotlib se renderizaran directamente en nuestro cuaderno. También hay muchas otras magias que podemos usar.

Cómo usar magia en Jupyter

Un buen primer paso es abrir un Jupyter Notebook, escribir %lsmagicen una celda y ejecutar la celda. Esto generará una lista de las magias de línea y de celda disponibles, y también te dirá si “automagic” está activado. 

  • Las magias de línea operan en una sola línea de una celda de código
  • La magia celular opera en toda la celda de código en la que se llaman

Si automagic está activado, puede ejecutar una magia simplemente escribiéndola en su propia línea en una celda de código y ejecutando la celda. Si está desactivado, deberá poner   %antes de las magias de línea y   %%  antes de las magias de células para usarlas.

Muchas magias requieren información adicional (al igual que una función requiere un argumento) para decirles cómo operar. Veremos un ejemplo en la siguiente sección, pero puede ver la documentación de cualquier magia ejecutándola con un signo de interrogación, así:

%matplotlib?

Cuando ejecute la celda anterior en un cuaderno, aparecerá una larga cadena de documentos en la pantalla con detalles sobre cómo puede usar la magia.

Algunos comandos mágicos útiles

Cubrimos más en el tutorial avanzado de Jupyter , pero aquí hay algunos para comenzar:

Mando mágicoQue hace
%correrEjecuta un archivo de secuencia de comandos externo como parte de la celda que se está ejecutando.
Por ejemplo, si aparece % run myscript.py en una celda de código, el kernel ejecutará myscript.py como parte de esa celda.
%cronométraloCuenta bucles, mide e informa cuánto tiempo tarda en ejecutarse una celda de código.
% archivo de escrituraGuarde el contenido de una celda en un archivo.
Por ejemplo,  % savefile myscript.py guardaría la celda de código como un archivo externo llamado myscript.py.
%TiendaGuarde una variable para usarla en un cuaderno diferente.
% pwdImprima la ruta del directorio en el que está trabajando actualmente.
%% javascriptEjecuta la celda como código JavaScript.

Hay mucho más de donde vino eso. Ingrese a Jupyter Notebooks y comience a explorar usando%lsmagic !

Pensamientos finales

Comenzando desde cero, nos hemos familiarizado con el flujo de trabajo natural de Jupyter Notebooks, profundizamos en las funciones más avanzadas de IPython y finalmente aprendimos cómo compartir nuestro trabajo con amigos, colegas y el mundo. ¡Y logramos todo esto desde un propio cuaderno!

Debe quedar claro cómo los cuadernos promueven una experiencia de trabajo productiva al reducir el cambio de contexto y emular un desarrollo natural de pensamientos durante un proyecto. El poder de usar Jupyter Notebooks también debería ser evidente, y cubrimos muchos clientes potenciales para que comience a explorar más funciones avanzadas en sus propios proyectos.

Si desea obtener más inspiración para sus propios cuadernos, Jupyter ha creado  una galería de cuadernos de Jupyter interesantes  que pueden resultarle útiles y la  página de inicio de Nbviewer.  enlaces a algunos ejemplos realmente elegantes de cuadernos de calidad.

Estimación de máxima verosimilitud en Python


¿Qué pasa si una relación lineal no es un supuesto apropiado para nuestro modelo?

Una alternativa ampliamente propuesta por Thomas J. Sargent y John Stachurski es la estimación de máxima verosimilitud, que implica especificar una clase de distribuciones, indexadas por parámetros desconocidos, y luego usar los datos para precisar los valores de estos parámetros.

El beneficio relativo a la regresión lineal es que permite una mayor flexibilidad en las relaciones probabilísticas entre variables.

Aquí ilustramos la máxima probabilidad replicando el artículo de Daniel Treisman (2016) sobre los multimillonarios de Rusia , que conecta la cantidad de multimillonarios en un país con sus características económicas concluyendo que Rusia tiene una mayor cantidad de multimillonarios de lo que predicen factores económicos como el tamaño del mercado y la tasa impositiva.

Requeriremos las siguientes importaciones:

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
import numpy as np
from numpy import exp
from scipy.special import factorial
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import statsmodels.api as sm
from statsmodels.api import Poisson
from scipy import stats
from scipy.stats import norm
from statsmodels.iolib.summary2 import summary_col

 Configuración y suposiciones 

Consideremos los pasos que debemos seguir en la estimación de máxima verosimilitud y cómo pertenecen a este estudio.

Flujo de ideas 

El primer paso con la estimación de máxima verosimilitud es elegir la distribución de probabilidad que se cree que genera los datos.

Más precisamente, necesitamos suponer qué clase de distribución paramétrica está generando los datos.

  • por ejemplo, la clase de todas las distribuciones normales o la clase de todas las distribuciones gamma.

Cada una de estas clases es una familia de distribuciones indexadas por un número finito de parámetros.

  • Por ejemplo, la clase de distribuciones normales es una familia de distribuciones indexadas por su media. μ∈(−∞,∞) y desviación estándar σ∈(0,∞).

Dejaremos que los datos seleccionen un elemento particular de la clase definiendo los parámetros.

Las estimaciones de parámetros así producidas se denominarán estimaciones de máxima verosimilitud .

 Contando multimillonarios 

Treisman está interesado en estimar el número de multimillonarios en diferentes países.El número de multimillonarios se valora en números enteros.Por lo tanto, consideramos distribuciones que toman valores solo en los enteros no negativos.

(Esta es una de las razones por las que la regresión por mínimos cuadrados no es la mejor herramienta para el problema actual, ya que la variable dependiente en la regresión lineal no está restringida a valores enteros)

Una distribución entera es la distribución de Poisson , cuya función de masa de probabilidad (pmf) esf(y)=μyy!e−μ,y=0,1,2,…,∞

Podemos trazar la distribución de Poisson sobre y para diferentes valores de μ como sigue

poisson_pmf = lambda y, μ: μ**y / factorial(y) * exp(-μ)
y_values = range(0, 25)

fig, ax = plt.subplots(figsize=(12, 8))

for μ in [1, 5, 10]:
    distribution = []
    for y_i in y_values:
        distribution.append(poisson_pmf(y_i, μ))
    ax.plot(y_values,
            distribution,
            label=f'$\mu$={μ}',
            alpha=0.5,
            marker='o',
            markersize=8)

ax.grid()
ax.set_xlabel('$y$', fontsize=14)
ax.set_ylabel('$f(y \mid \mu)$', fontsize=14)
ax.axis(xmin=0, ymin=0)
ax.legend(fontsize=14)

plt.show()

_images / mle_3_0.png

Observe que la distribución de Poisson comienza a parecerse a una distribución normal como la media de y aumenta.

Echemos un vistazo a la distribución de los datos con los que trabajaremos en esta conferencia.

La principal fuente de datos de Treisman son las clasificaciones anuales de multimillonarios de Forbes y su patrimonio neto estimado.

El conjunto de datos mle/fp.dtase puede descargar desde aquí o desde su página AER .

pd.options.display.max_columns = 10

# Load in data and view
df = pd.read_stata('https://github.com/QuantEcon/lecture-python/blob/master/source/_static/lecture_specific/mle/fp.dta?raw=true')
df.head()

paísccodeañoañoentumecidotopint08rintrnoyrsroflawnrrents
0Estados Unidos2.01990.021990.0Yaya39.7999994.98840520,01,61Yaya
1Estados Unidos2.01991.021991.0Yaya39.7999994.98840520,01,61Yaya
2Estados Unidos2.01992,021992.0Yaya39.7999994.98840520,01,61Yaya
3Estados Unidos2.01993.021993.0Yaya39.7999994.98840520,01,61Yaya
4Estados Unidos2.01994.021994.0Yaya39.7999994.98840520,01,61Yaya

5 filas × 36 columnas

Usando un histograma, podemos ver la distribución del número de multimillonarios por país numbil0, en 2008 (los Estados Unidos se descartan para fines de trazado)

numbil0_2008 = df[(df['year'] == 2008) & (
    df['country'] != 'United States')].loc[:, 'numbil0']

plt.subplots(figsize=(12, 8))
plt.hist(numbil0_2008, bins=30)
plt.xlim(left=0)
plt.grid()
plt.xlabel('Number of billionaires in 2008')
plt.ylabel('Count')
plt.show()

_images / mle_7_0.png

A partir del histograma, parece que la suposición de Poisson no es irrazonable (aunque con un valor muy bajo μ y algunos valores atípicos).

Distribuciones condicionales 

En el artículo de Treisman, la variable dependiente: el número de multimillonarios yi en el pais i – se modela en función del PIB per cápita, el tamaño de la población y los años de pertenencia al GATT y la OMC.

Por tanto, la distribución de yi necesita estar condicionado al vector de variables explicativas xi.

La formulación estándar, el llamado modelo de regresión de Poisson , es la siguiente:(58.1)f(yi∣xi)=μiyiyi!e−μi;yi=0,1,2,…,∞.where μi=exp⁡(xi′β)=exp⁡(β0+β1xi1+…+βkxik)

Para ilustrar la idea de que la distribución de yi depende de xi ejecutemos una simulación simple.

Usamos nuestra poisson_pmffunción de arriba y valores arbitrarios para β y xi

y_values = range(0, 20)

# Define a parameter vector with estimates
β = np.array([0.26, 0.18, 0.25, -0.1, -0.22])

# Create some observations X
datasets = [np.array([0, 1, 1, 1, 2]),
            np.array([2, 3, 2, 4, 0]),
            np.array([3, 4, 5, 3, 2]),
            np.array([6, 5, 4, 4, 7])]


fig, ax = plt.subplots(figsize=(12, 8))

for X in datasets:
    μ = exp(X @ β)
    distribution = []
    for y_i in y_values:
        distribution.append(poisson_pmf(y_i, μ))
    ax.plot(y_values,
            distribution,
            label=f'$\mu_i$={μ:.1}',
            marker='o',
            markersize=8,
            alpha=0.5)

ax.grid()
ax.legend()
ax.set_xlabel('$y \mid x_i$')
ax.set_ylabel(r'$f(y \mid x_i; \beta )$')
ax.axis(xmin=0, ymin=0)
plt.show()

_images / mle_9_0.png

Podemos ver que la distribución de yi está condicionado a xi (μi ya no es constante).

Estimación de máxima verosimilitud 

En nuestro modelo para el número de multimillonarios, la distribución condicional contiene 4 (k=4) parámetros que necesitamos estimar.

Etiquetaremos todo nuestro vector de parámetros como β dóndeβ=[β0β1β2β3]

Para estimar el modelo usando MLE, queremos maximizar la probabilidad de que nuestra estimación β^ es el verdadero parámetro β.

Intuitivamente, queremos encontrar el β^ que mejor se adapte a nuestros datos.

Primero, necesitamos construir la función de verosimilitud L(β), que es similar a una función de densidad de probabilidad conjunta.

Supongamos que tenemos algunos datos yi={y1,y2} y yi∼f(yi).

Si y1 y y2 son independientes, la PMF conjunta de estos datos es f(y1,y2)=f(y1)⋅f(y2).

Si yi sigue una distribución de Poisson con λ=7, podemos visualizar el PMF conjunto así

def plot_joint_poisson(μ=7, y_n=20):
    yi_values = np.arange(0, y_n, 1)

    # Create coordinate points of X and Y
    X, Y = np.meshgrid(yi_values, yi_values)

    # Multiply distributions together
    Z = poisson_pmf(X, μ) * poisson_pmf(Y, μ)

    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(X, Y, Z.T, cmap='terrain', alpha=0.6)
    ax.scatter(X, Y, Z.T, color='black', alpha=0.5, linewidths=1)
    ax.set(xlabel='$y_1$', ylabel='$y_2$')
    ax.set_zlabel('$f(y_1, y_2)$', labelpad=10)
    plt.show()

plot_joint_poisson(μ=7, y_n=20)

_images / mle_11_0.png

De manera similar, la pmf conjunta de nuestros datos (que se distribuye como una distribución condicional de Poisson) se puede escribir comof(y1,y2,…,yn∣x1,x2,…,xn;β)=∏i=1nμiyiyi!e−μi

yi está condicionado a ambos valores de xi y los parámetros β.

La función de verosimilitud es la misma que la pmf conjunta, pero trata el parámetro β como una variable aleatoria y toma las observaciones (yi,xi) como se indicaL(β∣y1,y2,…,yn ; x1,x2,…,xn)=∏i=1nμiyiyi!e−μi=f(y1,y2,…,yn∣ x1,x2,…,xn;β)

Ahora que tenemos nuestra función de verosimilitud, queremos encontrar el β^ que produce el valor de máxima verosimilitudmaxβL(β)

Al hacerlo, generalmente es más fácil maximizar la probabilidad logarítmica (considere diferenciar f(x)=xexp⁡(x) vs. f(x)=log⁡(x)+x).

Dado que tomar un logaritmo es una transformación creciente monótona, un maximizador de la función de verosimilitud también será un maximizador de la función logarítmica de verosimilitud.

En nuestro caso, la probabilidad logarítmica eslog⁡L(β)= log⁡(f(y1;β)⋅f(y2;β)⋅…⋅f(yn;β))=∑i=1nlog⁡f(yi;β)=∑i=1nlog⁡(μiyiyi!e−μi)=∑i=1nyilog⁡μi−∑i=1nμi−∑i=1nlog⁡y!

El MLE del Poisson al Poisson para β^ se puede obtener resolviendomaxβ(∑i=1nyilog⁡μi−∑i=1nμi−∑i=1nlog⁡y!)

Sin embargo, no existe una solución analítica para el problema anterior: para encontrar el MLE necesitamos usar métodos numéricos.

 MLE con métodos numéricos 

Muchas distribuciones no tienen buenas soluciones analíticas y, por lo tanto, requieren métodos numéricos para resolver las estimaciones de los parámetros.

Uno de esos métodos numéricos es el algoritmo de Newton-Raphson.

Nuestro objetivo es encontrar la estimación de máxima verosimilitud β^.

A β^, la primera derivada de la función logarítmica de verosimilitud será igual a 0.

Ilustremos esto suponiendolog⁡L(β)=−(β−10)2−10

β = np.linspace(1, 20)
logL = -(β - 10) ** 2 - 10
dlogL = -2 * β + 20

fig, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(12, 8))

ax1.plot(β, logL, lw=2)
ax2.plot(β, dlogL, lw=2)

ax1.set_ylabel(r'$log \mathcal{L(\beta)}$',
               rotation=0,
               labelpad=35,
               fontsize=15)
ax2.set_ylabel(r'$\frac{dlog \mathcal{L(\beta)}}{d \beta}$ ',
               rotation=0,
               labelpad=35,
               fontsize=19)
ax2.set_xlabel(r'$\beta$', fontsize=15)
ax1.grid(), ax2.grid()
plt.axhline(c='black')
plt.show()

_images / mle_13_0.png

La gráfica muestra que el valor de máxima verosimilitud (la gráfica superior) ocurre cuando dlog⁡L(β)dβ=0 (la trama inferior).

Por lo tanto, la probabilidad se maximiza cuando β=10.

También podemos asegurarnos de que este valor sea un máximo (en lugar de un mínimo) comprobando que la segunda derivada (pendiente del gráfico inferior) sea negativa.

El algoritmo de Newton-Raphson encuentra un punto donde la primera derivada es 0.

Para usar el algoritmo, hacemos una estimación inicial del valor máximo, β0 (las estimaciones de los parámetros de MCO pueden ser una suposición razonable), entonces

  1. Utilice la regla de actualización para iterar el algoritmoβ(k+1)=β(k)−H−1(β(k))G(β(k))dónde:G(β(k))=dlog⁡L(β(k))dβ(k)H(β(k))=d2log⁡L(β(k))dβ(k)dβ(k)′
  2. Compruebe si β(k+1)−β(k)<tol
    • Si es verdadero, deje de iterar y configure β^=β(k+1)
    • Si es falso, actualice β(k+1)

Como puede verse en la ecuación de actualización, β(k+1)=β(k) sólo cuando G(β(k))=0es decir. donde la primera derivada es igual a 0.

(En la práctica, dejamos de iterar cuando la diferencia está por debajo de un pequeño umbral de tolerancia)

Intentemos implementar el algoritmo de Newton-Raphson.

Primero, crearemos una clase llamada PoissonRegressionpara que podamos volver a calcular fácilmente los valores de probabilidad de registro, gradiente y hessiano para cada iteración.

class PoissonRegression:

    def __init__(self, y, X, β):
        self.X = X
        self.n, self.k = X.shape
        # Reshape y as a n_by_1 column vector
        self.y = y.reshape(self.n,1)
        # Reshape β as a k_by_1 column vector
        self.β = β.reshape(self.k,1)

    def μ(self):
        return np.exp(self.X @ self.β)

    def logL(self):
        y = self.y
        μ = self.μ()
        return np.sum(y * np.log(μ) - μ - np.log(factorial(y)))

    def G(self):
        y = self.y
        μ = self.μ()
        return X.T @ (y - μ)

    def H(self):
        X = self.X
        μ = self.μ()
        return -(X.T @ (μ * X))

Nuestra función newton_raphsontomará un PoissonRegressionobjeto que tiene una suposición inicial del vector de parámetrosβ0.

El algoritmo actualizará el vector de parámetros de acuerdo con la regla de actualización y volverá a calcular el gradiente y las matrices hessianas en las estimaciones de los nuevos parámetros.

La iteración terminará cuando:

  • La diferencia entre el parámetro y el parámetro actualizado está por debajo de un nivel de tolerancia.
  • Se ha alcanzado el número máximo de iteraciones (lo que significa que no se logra la convergencia).

Para que podamos tener una idea de lo que sucede mientras se ejecuta el algoritmo, display=Truese agrega una opción para imprimir valores en cada iteración.

def newton_raphson(model, tol=1e-3, max_iter=1000, display=True):

    i = 0
    error = 100  # Initial error value

    # Print header of output
    if display:
        header = f'{"Iteration_k":<13}{"Log-likelihood":<16}{"θ":<60}'
        print(header)
        print("-" * len(header))

    # While loop runs while any value in error is greater
    # than the tolerance until max iterations are reached
    while np.any(error > tol) and i < max_iter:
        H, G = model.H(), model.G()
        β_new = model.β - (np.linalg.inv(H) @ G)
        error = β_new - model.β
        model.β = β_new

        # Print iterations
        if display:
            β_list = [f'{t:.3}' for t in list(model.β.flatten())]
            update = f'{i:<13}{model.logL():<16.8}{β_list}'
            print(update)

        i += 1

    print(f'Number of iterations: {i}')
    print(f'β_hat = {model.β.flatten()}')

    # Return a flat array for β (instead of a k_by_1 column vector)
    return model.β.flatten()

Probemos nuestro algoritmo con un pequeño conjunto de datos de 5 observaciones y 3 variables en X.

X = np.array([[1, 2, 5],
              [1, 1, 3],
              [1, 4, 2],
              [1, 5, 2],
              [1, 3, 1]])

y = np.array([1, 0, 1, 1, 0])

# Take a guess at initial βs
init_β = np.array([0.1, 0.1, 0.1])

# Create an object with Poisson model values
poi = PoissonRegression(y, X, β=init_β)

# Use newton_raphson to find the MLE
β_hat = newton_raphson(poi, display=True)

Iteration_k  Log-likelihood  θ                                                           
-----------------------------------------------------------------------------------------
0            -4.3447622      ['-1.49', '0.265', '0.244']
1            -3.5742413      ['-3.38', '0.528', '0.474']
2            -3.3999526      ['-5.06', '0.782', '0.702']
3            -3.3788646      ['-5.92', '0.909', '0.82']
4            -3.3783559      ['-6.07', '0.933', '0.843']
5            -3.3783555      ['-6.08', '0.933', '0.843']
Number of iterations: 6
β_hat = [-6.07848205  0.93340226  0.84329625]

Como se trataba de un modelo simple con pocas observaciones, el algoritmo logró la convergencia en solo 6 iteraciones.

Puede ver que con cada iteración, el valor de probabilidad logarítmica aumentó.

Recuerde, nuestro objetivo era maximizar la función de probabilidad logarítmica, que el algoritmo ha trabajado para lograr.

Además, tenga en cuenta que el aumento en log⁡L(β(k)) se vuelve más pequeño con cada iteración.

Esto se debe a que el gradiente se acerca a 0 a medida que alcanzamos el máximo y, por lo tanto, el numerador de nuestra ecuación de actualización se hace más pequeño.

El vector de gradiente debe estar cerca de 0 en β^

poi.G()

array([[-3.95169228e-07],
       [-1.00114805e-06],
       [-7.73114562e-07]])

El proceso iterativo se puede visualizar en el siguiente diagrama, donde el máximo se encuentra en β=10

logL = lambda x: -(x - 10) ** 2 - 10

def find_tangent(β, a=0.01):
    y1 = logL(β)
    y2 = logL(β+a)
    x = np.array([[β, 1], [β+a, 1]])
    m, c = np.linalg.lstsq(x, np.array([y1, y2]), rcond=None)[0]
    return m, c

β = np.linspace(2, 18)
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(β, logL(β), lw=2, c='black')

for β in [7, 8.5, 9.5, 10]:
    β_line = np.linspace(β-2, β+2)
    m, c = find_tangent(β)
    y = m * β_line + c
    ax.plot(β_line, y, '-', c='purple', alpha=0.8)
    ax.text(β+2.05, y[-1], f'$G({β}) = {abs(m):.0f}$', fontsize=12)
    ax.vlines(β, -24, logL(β), linestyles='--', alpha=0.5)
    ax.hlines(logL(β), 6, β, linestyles='--', alpha=0.5)

ax.set(ylim=(-24, -4), xlim=(6, 13))
ax.set_xlabel(r'$\beta$', fontsize=15)
ax.set_ylabel(r'$log \mathcal{L(\beta)}$',
               rotation=0,
               labelpad=25,
               fontsize=15)
ax.grid(alpha=0.3)
plt.show()

_images / mle_23_0.png

Tenga en cuenta que nuestra implementación del algoritmo de Newton-Raphson es bastante básica; para implementaciones más sólidas, consulte, por ejemplo, scipy.optimize .

Estimación de máxima verosimilitud con statsmodels

Ahora que sabemos lo que está pasando bajo el capó, podemos aplicar MLE a una aplicación interesante.

Usaremos el modelo de regresión de Poisson statsmodelspara obtener una salida más rica con errores estándar, valores de prueba y más.

statsmodels utiliza el mismo algoritmo anterior para encontrar las estimaciones de máxima verosimilitud.

Antes de comenzar, volvamos a estimar nuestro modelo simple con statsmodels para confirmar que obtenemos los mismos coeficientes y valor logarítmico de verosimilitud.

X = np.array([[1, 2, 5],
              [1, 1, 3],
              [1, 4, 2],
              [1, 5, 2],
              [1, 3, 1]])

y = np.array([1, 0, 1, 1, 0])

stats_poisson = Poisson(y, X).fit()
print(stats_poisson.summary())

Optimization terminated successfully.
         Current function value: 0.675671
         Iterations 7
                          Poisson Regression Results                          
==============================================================================
Dep. Variable:                      y   No. Observations:                    5
Model:                        Poisson   Df Residuals:                        2
Method:                           MLE   Df Model:                            2
Date:                Thu, 15 Jul 2021   Pseudo R-squ.:                  0.2546
Time:                        00:14:01   Log-Likelihood:                -3.3784
converged:                       True   LL-Null:                       -4.5325
Covariance Type:            nonrobust   LLR p-value:                    0.3153
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -6.0785      5.279     -1.151      0.250     -16.425       4.268
x1             0.9334      0.829      1.126      0.260      -0.691       2.558
x2             0.8433      0.798      1.057      0.291      -0.720       2.407
==============================================================================

Ahora repliquemos los resultados del artículo de Daniel Treisman, Los multimillonarios de Rusia , mencionado anteriormente en la conferencia.

Treisman comienza estimando la ecuación (58.1) , donde:

  • yi es number of billionairesi
  • xi1 es log⁡GDP per capitai
  • xi2 es log⁡populationi
  • xi3 es years in GATTi – años de membresía en el GATT y la OMC (para poder acceder a los mercados internacionales)

El documento solo considera el año 2008 para la estimación.

Configuraremos nuestras variables para la estimación de esa manera (debe tener los datos asignados dfdesde antes en la lección)

# Keep only year 2008
df = df[df['year'] == 2008]

# Add a constant
df['const'] = 1

# Variable sets
reg1 = ['const', 'lngdppc', 'lnpop', 'gattwto08']
reg2 = ['const', 'lngdppc', 'lnpop',
        'gattwto08', 'lnmcap08', 'rintr', 'topint08']
reg3 = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08',
        'rintr', 'topint08', 'nrrents', 'roflaw']

Entonces podemos usar la Poissonfunción de statsmodelspara ajustar el modelo.

Usaremos errores estándar robustos como en el artículo del autor.

# Specify model
poisson_reg = sm.Poisson(df[['numbil0']], df[reg1],
                         missing='drop').fit(cov_type='HC0')
print(poisson_reg.summary())

Optimization terminated successfully.
         Current function value: 2.226090
         Iterations 9
                          Poisson Regression Results                          
==============================================================================
Dep. Variable:                numbil0   No. Observations:                  197
Model:                        Poisson   Df Residuals:                      193
Method:                           MLE   Df Model:                            3
Date:                Thu, 15 Jul 2021   Pseudo R-squ.:                  0.8574
Time:                        00:14:01   Log-Likelihood:                -438.54
converged:                       True   LL-Null:                       -3074.7
Covariance Type:                  HC0   LLR p-value:                     0.000
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const        -29.0495      2.578    -11.268      0.000     -34.103     -23.997
lngdppc        1.0839      0.138      7.834      0.000       0.813       1.355
lnpop          1.1714      0.097     12.024      0.000       0.980       1.362
gattwto08      0.0060      0.007      0.868      0.386      -0.008       0.019
==============================================================================

¡Éxito! El algoritmo pudo lograr la convergencia en 9 iteraciones.

Nuestro resultado indica que el PIB per cápita, la población y los años de membresía en el Acuerdo General sobre Aranceles Aduaneros y Comercio (GATT) están relacionados positivamente con la cantidad de multimillonarios que tiene un país, como se esperaba.

Estimemos también los modelos más completos del autor y los mostraremos en una sola tabla.

regs = [reg1, reg2, reg3]
reg_names = ['Model 1', 'Model 2', 'Model 3']
info_dict = {'Pseudo R-squared': lambda x: f"{x.prsquared:.2f}",
             'No. observations': lambda x: f"{int(x.nobs):d}"}
regressor_order = ['const',
                   'lngdppc',
                   'lnpop',
                   'gattwto08',
                   'lnmcap08',
                   'rintr',
                   'topint08',
                   'nrrents',
                   'roflaw']
results = []

for reg in regs:
    result = sm.Poisson(df[['numbil0']], df[reg],
                        missing='drop').fit(cov_type='HC0',
                                            maxiter=100, disp=0)
    results.append(result)

results_table = summary_col(results=results,
                            float_format='%0.3f',
                            stars=True,
                            model_names=reg_names,
                            info_dict=info_dict,
                            regressor_order=regressor_order)
results_table.add_title('Table 1 - Explaining the Number of Billionaires \
                        in 2008')
print(results_table)

Table 1 - Explaining the Number of Billionaires                         in 2008
=================================================
                  Model 1    Model 2    Model 3  
-------------------------------------------------
const            -29.050*** -19.444*** -20.858***
                 (2.578)    (4.820)    (4.255)   
lngdppc          1.084***   0.717***   0.737***  
                 (0.138)    (0.244)    (0.233)   
lnpop            1.171***   0.806***   0.929***  
                 (0.097)    (0.213)    (0.195)   
gattwto08        0.006      0.007      0.004     
                 (0.007)    (0.006)    (0.006)   
lnmcap08                    0.399**    0.286*    
                            (0.172)    (0.167)   
rintr                       -0.010     -0.009    
                            (0.010)    (0.010)   
topint08                    -0.051***  -0.058*** 
                            (0.011)    (0.012)   
nrrents                                -0.005    
                                       (0.010)   
roflaw                                 0.203     
                                       (0.372)   
Pseudo R-squared 0.86       0.90       0.90      
No. observations 197        131        131       
=================================================
Standard errors in parentheses.
* p<.1, ** p<.05, ***p<.01

El resultado sugiere que la frecuencia de multimillonarios se correlaciona positivamente con el PIB per cápita, el tamaño de la población, la capitalización del mercado de valores y se correlaciona negativamente con la tasa de impuesto sobre la renta marginal máxima.

Para analizar nuestros resultados por país, podemos trazar la diferencia entre los valores predichos y reales, luego ordenar de mayor a menor y trazar los primeros 15

data = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08', 'rintr',
        'topint08', 'nrrents', 'roflaw', 'numbil0', 'country']
results_df = df[data].dropna()

# Use last model (model 3)
results_df['prediction'] = results[-1].predict()

# Calculate difference
results_df['difference'] = results_df['numbil0'] - results_df['prediction']

# Sort in descending order
results_df.sort_values('difference', ascending=False, inplace=True)

# Plot the first 15 data points
results_df[:15].plot('country', 'difference', kind='bar',
                    figsize=(12,8), legend=False)
plt.ylabel('Number of billionaires above predicted level')
plt.xlabel('Country')
plt.show()

_images / mle_33_0.png

Como podemos ver, Rusia tiene, con mucho, el mayor número de multimillonarios por encima de lo que predice el modelo (alrededor de 50 más de lo esperado).

Treisman usa este resultado empírico para discutir las posibles razones del exceso de multimillonarios de Rusia, incluido el origen de la riqueza en Rusia, el clima político y la historia de la privatización en los años posteriores a la URSS.

Fuente https://python.quantecon.org/mle.html

Introducción a python científico


Muchas áreas de carácter científico-técnico la adecuada elección del software y/o lenguaje de programación empleado es determinante, de cara a la potencia, versatilidad, facilidad de uso y acceso por parte de todos los usuarios en sus propios dispositivos, de manera generalizada y gratuita.

Dentro del software libre, uno de los que últimamente ha tenido una mejora sustancial, con la inclusión de potentes y versátiles nuevos módulos de cálculo simbólico (SymPy), numérico (NumPy, SciPy) y gráfico (PyPlot y Matplotlib) ha sido sin duda Python, y de ahí su vertiginosa evolución y expansión a nivel mundial, no sólo en el ámbito académico, sino también en el científico e industrial. De hecho, basta con echar un vistazo a las numerosas propuestas, tanto de comunidades de desarrolladores como de empresas privadas, surgidas a raíz de la versión de base inicial de Python, como por ejemplo IPython (interface interactivo de fácil uso, que gracias a Jupyter Notebook permite una versión HTML similar a los notebooks de Mathematica o Mapple) o Spyder (entorno integrado para cálculo científico parecido al de Matlab u Octave).

Por otro lado existen versiones completas de desarrollo, integrando Python como soporte de cálculo, pero con editores avanzados de texto para la programación y la depuración de código, ventanas de gráficos y datos, etc. La mayoría de estas plataformas integradas están disponibles para los distintos sistemas operativos Linux, MacOS X y Windows. Entre ellas cabría destacar Enthought Python Distribution (EPD), PyCharm y Anaconda CE (de Continuum Analytics).

Aunque no podamos abarcar todos los aspectos “básicos” del python científico, intentaremos en este resumen dar una idea de las principales librerías un funciones que podemos usar para NILM (Non-Intrusive Load Monitoring) sin olvidar los fundamentos de :Matplotlib y Numpy

Matplotlib: visualización con Python 

Matplotlib es una biblioteca completa para crear visualizaciones estáticas, animadas e interactivas en Python haciendo que las cosas fáciles sean fáciles y las difíciles posibles.

Nos permite crear :

  • Desarrollando gráficos de calidad de publicación con solo unas pocas líneas de código
  • Utilizando figuras interactivas que puedan hacer zoom, desplazarse, actualizar …

Personalizar

  • Tomando el control total de los estilos de línea, las propiedades de la fuente, las propiedades de los ejes …
  • Exportando e incrustando en varios formatos de archivo y entornos interactivos

Ampliar

  • Explorando la funcionalidad personalizada proporcionada por paquetes de terceros
  • Obteniendo más información sobre Matplotlib a través de los numerosos recursos de aprendizaje externos

Matplotlib es en resumen la librería de python para dibujar (equivalente al plot en matlab).

matplotlib

Puede encontrar mas información en el sitio oficial https://matplotlib.org/

Numpy

NumPy es una biblioteca para el lenguaje de programación Python que da soporte para crear vectores y matrices grandes multidimensionales, junto con una gran colección de funciones matemáticas de alto nivel para operar con ellas.

Numpy es pues una librería especializada para operaciones con matrices y vectores

Puede encontrar mas información en l sitio oficial https://numpy.org/

La imagen tiene un atributo ALT vacío; su nombre de archivo es image-11.png

Primeros pasos

Primero, es necesario importarlas al workspace

import numpy as np
import matplotlib.pyplot as plt

Opciones de visualizacion de matplotlib para un notebook

%matplotlib inline
plt.rcParams['figure.figsize'] = (13, 6)
plt.style.use('ggplot')

Otras importaciones:

import warnings
warnings.filterwarnings('ignore')

Crear arrays en python es muy sencillo y se puede hacer de forma nativa usando un tipo list. Sin embargo, aquí consideramos arrays del tipo numpy pues esto arrays incluyen funciones que facilitan las operaciones matemáticas y su manipulación

v=[1,2,3] # tipo list
v=np.array([1,2,3]) # array numpy
print (v)
print ("Dimensiones: " + str(v.ndim)) # numero de dimensiones
print ("Elementos: " + str(v.size)) # numero de elementos
print ("Longitud de las dimensiones: " + str(v.shape)) # longitud de cada dimensión
[1 2 3]
Dimensiones: 1
Elementos: 3
Longitud de las dimensiones: (3,)

Crear una matriz de 2 x 3:

v=np.array([[1,2,3], [4,5,6]])
print (v)
print ('Dimensiones: ' + str(v.ndim)) # numero de dimensiones
print ('Elementos: '+str(v.size)) # numero de elementos
print ('Longitud de las dimensiones: '+str(v.shape)) # longitud de cada dimensión
[[1 2 3]
 [4 5 6]]
Dimensiones: 2
Elementos: 6
Longitud de las dimensiones: (2, 3)

Crear una Matriz triple de 2 x 3 x 2 :

v=np.array([[[1,2], [3,4]],[[5,6],  [7,8]]])
print (v)
print ("Dimensiones: " + str(v.ndim)) # numero de dimensiones
print ("Elementos: "+str(v.size)) # numero de elementos
print ("Longitud de las dimensiones: "+str(v.shape) )# longitud de cada dimensión

[[[1 2], 
      [3 4]],

 [[5 6], 
      [7 8]]]
Dimensiones: 3
Elementos: 8
Longitud de las dimensiones: (2, 2, 2)

Utilizamos la función reshape para redimensionar los arrays

1 dimension

print (v.reshape(8,))
[1 2 3 4 5 6 7 8]

2 dimensiones

print (v.reshape(2,4))
[[1 2 3 4]
 [5 6 7 8]]

Matriz Identidad de 5×5

print (np.identity(5))
[[1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]

Matriz de unos de 5×5

print ( np.ones([5,5]))
[[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]

Matriz de ceros de 5×5:

print (np.zeros([5,5]))
[[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]

Las operaciones por definición son elementwise

a=np.arange(5)
b=10*np.ones(5)
print ("vector a: "+str(a))
print ("vector b: "+str(b))
print ("suma b+a: "+str(b-a))
print ("resta b-a: "+str(b+a))
print ("producto b*a: "+str(b*a))
vector a: [0 1 2 3 4]
vector b: [10. 10. 10. 10. 10.]
suma b+a: [10.  9.  8.  7.  6.]
resta b-a: [10. 11. 12. 13. 14.]
producto b*a: [ 0. 10. 20. 30. 40.]

El producto de los vectores es:

a.dot(b)
100.0

Para las matrices tenemos que:

a=np.identity(3)
b=np.array([[1,2,3],[4,5,6],[7,8,9]])
print ("matriz a:\n"+str(a))
print ("matriz b:\n"+str(b))
print ("producto a*b:\n"+str(a.dot(b)))
print ("producto elementwise a.*b:\n"+str(a*b))
matriz a:
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
matriz b:
[[1 2 3]
 [4 5 6]
 [7 8 9]]
producto a*b:
[[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
producto elementwise a.*b:
[[1. 0. 0.]
 [0. 5. 0.]
 [0. 0. 9.]]

Vector formado por un rango de valores:

print ("De 0 a 10: " + str(np.arange(10)))
print ("De 10 a 20 de paso 0.5: "+str(np.arange(10,20,0.5)))
De 0 a 10: [0 1 2 3 4 5 6 7 8 9]
De 10 a 20 de paso 0.5: [10.  10.5 11.  11.5 12.  12.5 13.  13.5 14.  14.5 15.  15.5 16.  16.5
 17.  17.5 18.  18.5 19.  19.5]

Función linspace:

np.linspace(0,2*np.pi,10) # de 0 a 2*pi en 10 puntos equidistantes
array([0.        , 0.6981317 , 1.3962634 , 2.0943951 , 2.7925268 ,
       3.4906585 , 4.1887902 , 4.88692191, 5.58505361, 6.28318531])

función random:

np.random.rand(10)
array([0.63623588, 0.83924558, 0.35833155, 0.33835148, 0.53247758,
       0.0950348 , 0.2805706 , 0.47285484, 0.8696919 , 0.78361161])

Dibujar una función seno

t = np.arange(0.0, 2.0, 0.01)
s = np.sin(2*np.pi*t)
plt.plot(t, s)
plt.xlabel('time (s)')
plt.ylabel('voltage (mV)')
plt.title('Sinusoidal')
plt.grid(True)

Dibujar una función chirp

x=np.linspace(0,3*np.pi,500)
plt.plot(x,np.sin(x**2))
plt.title("A simple chirp")
Text(0.5, 1.0, 'A simple chirp')

Definir una función en Python

En python, las funciones pueden estar definidas en cualquier parte pero siempre antes de su llamada. En python, las anidaciones (bucles, if conditions, functions, etc.) se realizan mediante indentación, no existe el statement end. Las funciones se definen así:

def funcion_suma(x): 
    suma=0
    for i in x: 
        suma=suma+i
    return suma 
v=np.arange(10)
print (funcion_suma(v))
45

Aunque, como hemos dicho antes, numpy facilita las operaciones matemáticas y ya incluye una serie de operaciones:

print (v.sum())
print (v.cumsum())
print (v.mean())
45
[ 0  1  3  6 10 15 21 28 36 45]
4.5

Para saber más sobre numpy:

https://docs.scipy.org/doc/numpy-dev/user/quickstart.html

http://www.sam.math.ethz.ch/~raoulb/teaching/PythonTutorial/intro_numpy.html

(O simplemente “googleando”: numpy tutorial)

Pandas

Pandas (Python Data Analysis Library) es una librería de python para el análisis y manipulación de una gran cantidad de datos. También facilita el uso de “timeseries”

La llamada a la librería es:

import pandas as pd

Dado un archivo csv, la función read_csv carga los datos en un dataframe

# El parámetro parse_dates indica a pandas que al cargar este csv la primera columna [0] es de tipo datetime
df=pd.read_csv('data/events.csv',parse_dates=[0])

Las primeras N filas del dataFrame se puede visualizar de la siguiente forma

N=4
df.head(N)
timestamplabelphase
02011-10-20 12:22:01.473111A
12011-10-20 12:37:40.507111A
22011-10-20 13:23:55.390111A
32011-10-20 13:39:08.157111A

Y las N últimas columnas

df.tail(N)
timestamplabelphase
24812011-10-27 12:57:17.079111A
24822011-10-27 13:10:45.112111A
24832011-10-27 13:54:08.862111A
24842011-10-27 14:07:21.612111A

Podemos filtar por un cierto valor

df[df.phase=='B'].head()
timestamplabelphase
122011-10-20 15:45:54.590204B
132011-10-20 15:47:31.223204B
142011-10-20 16:09:00.424204B
182011-10-20 17:42:00.657155B
192011-10-20 17:42:04.407157B

Y hacer agrupaciones

df2=df.sort_values(['label','phase']).groupby(['label','phase']).count()
df2
timestamp
labelphase
101B26
102B25
103B24
108A16
111A619
1125A1
1126A1
1127A1
1200A1
1201A1

161 rows × 1 columns

Las nuevas columnas se crean fácilmente. Compatible con numpy.

df['x']=25*np.random.rand(len(df))
df['y']=100*np.sin(2*np.pi*np.linspace(0,2*np.pi,len(df)))
df['z']=df.x+df.y
df.head(5)
timestamplabelphasexyz
02011-10-20 12:22:01.473111A0.0167760.0000000.016776
12011-10-20 12:37:40.507111A5.6400571.5892417.229298
22011-10-20 13:23:55.390111A5.6322523.1780818.810333
32011-10-20 13:39:08.157111A9.2871414.76611914.053259
42011-10-20 14:25:51.473111A9.3135696.35295215.666521

Para dibujar sólo necesitamos usar la función plot

df.z.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e3ba88080>

Existen ciertas funciones predefinidas que facilitan los cálculos

df.z.cumsum().plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e3b801e80>

Y se pueden concatenar

# Si integramos y derivamos obtenemos la misma señal
df.z.cumsum().diff().plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e32625320>

Unas de las herramientas más potentes de pandas es la manipulación de timeseries

df.index=df.timestamp
df.z.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e325e6518>

Podemos filtar por fechas

d1='2011-10-21'
d2='2011-10-23'
df[(df.index>d1)&(df.index<d2)].z.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e325d5588>

Existe una gran flexibilidad a la hora de resamplear un dataframe

# Cada día 
df.resample('1D',how='sum')
labelxyz
timestamp
2011-10-20233182044.60316112047.56801914092.171180
2011-10-21507894364.168081-4531.116500-166.948419
2011-10-221356456971.7315141811.9787858783.710300
2011-10-231021845445.726210-4029.5292421416.196968
2011-10-24462593554.5631237284.89716310839.460286
2011-10-25525523933.715959-4712.890577-779.174618
2011-10-26590843503.958137-7682.254604-4178.296467
2011-10-27181921325.9193697454.6146868780.534055
# Cada 6 horas
df.resample('6H',how='count')
timestamplabelphasexyz
timestamp
2011-10-20 12:00:00373737373737
2011-10-20 18:00:00135135135135135135
2011-10-21 00:00:00160160160160160160
2011-10-21 06:00:00676767676767
2011-10-21 12:00:00262626262626
2011-10-21 18:00:00828282828282
2011-10-22 00:00:00242424242424
2011-10-22 06:00:00107107107107107107
2011-10-22 12:00:00215215215215215215
2011-10-22 18:00:00203203203203203203
2011-10-23 00:00:00636363636363
2011-10-23 06:00:00646464646464
2011-10-23 12:00:00737373737373
2011-10-23 18:00:00237237237237237237
2011-10-24 00:00:00363636363636
2011-10-24 06:00:00393939393939
2011-10-24 12:00:00494949494949
2011-10-24 18:00:00162162162162162162
2011-10-25 00:00:00333333333333
2011-10-25 06:00:00919191919191
2011-10-25 12:00:00373737373737
2011-10-25 18:00:00152152152152152152
2011-10-26 00:00:00232323232323
2011-10-26 06:00:00616161616161
2011-10-26 12:00:00161616161616
2011-10-26 18:00:00196196196196196196
2011-10-27 00:00:00363636363636
2011-10-27 06:00:00555555555555
2011-10-27 12:00:00666666

Para aprender más sobre pandas:

http://pandas.pydata.org/pandas-docs/stable/tutorials.html

http://pandas.pydata.org/pandas-docs/stable/10min.html

Detector de eventos

Vamos a crear un detector de eventos.

Dado el consumo eléctrico de una vivienda (voltage y corriente) queremos detectar en que momento se produce una conexión de un dispositivo. Para ello, filtraremos la señal sinusoidal obteniendo el valor eficaz de la corriente cada cierto intervalo. Los cambios en el valor eficaz van a determinar las conexiones y desconexiones de los distintos dispositivos. Derivando este valor eficaz, obtenemos picos en los que existe un cambio en el valor eficaz y, por lo tanto, posibles candidatos a eventos de conexión/desconexión. Finalmente, usando un detector de picos filtraremos los eventos reales del resto.

Mediremos nuestros resultados usando métricas estándar de NILM.

Paso por paso

Importar pandas, numpy y matplotlib tal y como se ha visto anteriormente

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

Definir una funcion llamada rms_function que devuelva un valor rms y que tenga como parámetro de entrada un vector de valores

# función rms 
def rms_function(x): 
    return np.sqrt(np.mean(np.square(x)))

Usar el siguiente path para cargar los datos en un dataframe df de pandas. Como parámetros: el índice es la columna 0 (index_col) y la fecha está en la columna 1 (parse_dates)

path='data/smart_meter_data.csv'
df= ...
path='data/smart_meter_data.csv'
df=pd.read_csv(path, parse_dates=[1],index_col=[0])

Mostrar las 5 primeras columnas del dataframe

df.head(5)
datetimeivlabelappl_namephase
02011-10-20 12:21:58.9730000.444955159.194375111RefrigeratorA
12011-10-20 12:21:58.9730830.402501160.677554111RefrigeratorA
22011-10-20 12:21:58.9731660.444955161.845163111RefrigeratorA
32011-10-20 12:21:58.9732491.102993163.107443111RefrigeratorA
42011-10-20 12:21:58.9733321.952074164.243495111RefrigeratorA

Imprimir mínimo y máximo de datetime y la diferencia de ambos

print (df.datetime.min())
print (df.datetime.max())
print (df.datetime.max()-df.datetime.min())
2011-10-20 12:21:58.973000
2011-10-20 12:23:03.713996
0 days 00:01:04.740996

Seleccionar datetime como índice del dataframe df

df.index=df.datetime

Periodo y frequencia de muestreo

# frecuencia
ts=df.datetime.diff().mean().total_seconds()
print (str(ts)+' seconds')
fs=1/ts
print ( str(fs)+' Hz')
8.3e-05 seconds
12048.192771084337 Hz

Dibujar Voltage (v) haciendo zoom en el intervalo de 100ms (6 periodos aproximadamente)

d1='2011-10-20 12:22:29.9'
d2='2011-10-20 12:22:30'
df[(df.index>d1)&(df.index<d2)].v.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e2b9f86d8>
df.i.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e2e274d30>

Resamplear mediante la función resample de pandas a 50ms (’50L’). La función rms_function se pasará como parámetro para obtener cada valor del resampleado. El resultado debe de guardarse en un dataframe nuevo llamado rms . Dibujar el resultado.

rms=pd.DataFrame(df.i.resample(....))
rms=pd.DataFrame(df.i.resample('50L',how=rms_function))
rms.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e2f8848d0>

Hacer la derivada del dataframe rms y guardar el resultado en rms_diff.

rms_diff=rms.diff()

Dibujar el resultado (rms_diff)

rms_diff.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e2f7f9fd0>

Guardar los valores de la columna “i” en una variable “y” en forma de array

y=rms_diff.i.values
Modifica los parámetros th_noise y dist de la función detect_peaks para obener los índices de los eventos y evaluar las métricas. Realizar el proceso 3 veces. ¿ Con qué valores de th_noise y dist se obtienen mejores resultados en las métricas?
th_noise=5
dist=5
from detect_peaks import detect_peaks
indexes=detect_peaks(y,mph=th_noise,mpd=dist)
dates=rms_diff.ix[indexes].index

Cuantos eventos hemos detectado

print (str(len(indexes))+' eventos detectados')
8 eventos detectados

Dibujamos los eventos y la corriente en una misma gráfica

plt.vlines(dates,-80,80)
df.i.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7f6e2f79ca58>

Métricas

from IPython.display import Image
Image(filename='metricas1.png')
Image(filename='metricas2.png')

Obtener las métricas: recall, precision y F1

FP=0.
P=9.
N=len(df)-P
TN=N-FP
P=9.
N=len(df)-P
TP=8.
FP=0.
FN=1.
TN=N-FP
recall=TP/(TP+FN)
precision=TP/(TP+FP)
F1=2*precision*recall/(precision+recall)
print (recall)
print (precision)
print (F1)
0.8888888888888888
1.0
0.9411764705882353

*Parámetros optimizados: * th_noise=0.1 y dist=5

*Con esto obtenemos : * recall=1, precision=1 y F1=1;