Estadística en Python: análisis de datos multidimensionales y regresión lineal (Parte IV)

Hasta ahora hemos tratado con una única variable por separado. Ahora vamos a ver qué podemos hacer con varias variables en la misma muestra (conjunto de datos). Nos interesará saber si están relacionadas o no (independencia o no). Si existe relación (estan correlacionadas) vamos a construir un modelo de regresión lineal.

Distribución conjunta de frecuencias

En el caso de dos variables, podemos construir una distribución conjunta de frecuencias. Se trata de una tabla de doble entrada donde cada dimensión corresponde a cada variable y el valor de las celdas representa la frecuencia del par. Para ello podemos usar crosstab también (de hecho, su uso original es este).

Ejemplo: En las votaciones a alcalde de la ciudad de Valladolid se presentaban Rafael, Silvia y Olga. Analiza los resultados e informa de quién fue el ganador de las elecciones. ¿Quién fue el candidato favorito en el barrio de La Rondilla?

Como podéis ver, un humando podría haber sacado estas conclusiones observando simplemente la tabla conjunta de frecuencias. ¿Quién tiene más votos en total? Rafael, con 6 en All (la suma de los distritos). ¿Quién ha sacado más votos en La Rondilla? Silvia, con 4 en la columna de La Rondilla. Por último, ¿votó más gente en el Centro o en La Rondilla? Votaron más en La Rondilla (8 votos), que en el Centro (7 votos).

A las frecuencias All se las llama comúnmente distribuciones marginales. Cuando discriminamos las frecuencias a un solo valor de una variable, se habla de distribuciones condicionadas, en este caso hemos usado la distribución de votos condicionada al distrito La Rondilla. Estas distribuciones son univariantes como habréis sospechado.

Gráfico XY o bivariante

Una manera muy útil de observar posibles correlaciones es con el gráfico XY, solamente disponible para distribuciones bivariantes. Cada observación se representa en el plano como un punto. En Matplotlib podemos dibujarlo con scatter.

Ejemplo: Represente el gráfico XY de las variables ingresos y gastos de las familias.

En la imagen podemos ver cada dato representado por un punto. En este ejemplo puede apreciarse como los puntos estan en torno a una línea recta invisible.

Covarianza

Para medir la relación entre dos variables podemos definir la covarianza:

\(
cov_{x,y}=\frac{\sum_{i=1}^{N}(x_{i}-\bar{x})(y_{i}-\bar{y})}{N}
\)

Pandas trae el método cov para calcular la matriz de covarianzas. De esta matriz, obtendremos el valor que nos interesa.

¿Y la covarianza qué nos dice? Por si mismo, bastante poco. Como mucho, si es positivo nos dice que se relacionarían de forma directa y si es negativa de forma inversa. Pero la covarianza está presente en muchas fórmulas.

Coeficiente de correlación lineal de Pearson

\(
r_{x,y}=\frac{cov_{x,y}}{s_{x}s_{y}}
\)

Uno de los métodos que usa la covarianza (aunque Pandas lo va a hacer solo) es el coeficiente de correlación lineal de Pearson. Cuanto más se acerque a 1 o -1 más correlacionadas están las variables. Su uso en Pandas es muy similar a la covarianza.

En este ejemplo concreto, el coeficiente de correlación de Pearson nos da 0.976175. Se trata de un valor lo suficientemente alto como para plantearnos una correlación lineal. Es decir, que pueda ser aproximado por una recta. Si este coeficiente es igual a 1 o -1, se puede decir que una variable es fruto de una transformación lineal de la otra.

Ajuste lineal

Vamos a intentar encontrar un modelo lineal que se ajuste a nuestras observaciones y nos permita hacer predicciones. Esta recta se llamará recta de regresión y se calcula de la siguiente forma:

\(
\hat{y}-\bar{y}=\frac{cov_{x,y}}{s_{x}^2}(x-\bar{x})
\)

Usando las funciones de varianza, media y covarianza Pandas no es muy complicado hacer una recta de regresión:

Que podemos probar visualmente:

Sin embargo SciPy ya nos trae un método que calcula la pendiente, la ordenada en el origen, el coeficiente de correlación lineal de Pearson y mucho más en un solo lugar. Es mucho más eficiente, se trata de linregress.

Además, para calcular los valores del gráfico, he usado vectorize de NumPy, que permite mapear los arrays nativos de NumPy. Más eficiente. Mismo resultado.

La ley de Ohm

¿Quién no ha oído hablar de la Ley de Ohm? Se trata de una ley que relaciona la diferencia de potencial con el amperaje dando lugar a la resistencia. La ley fue enunciada por George Simon Ohm, aunque no exactamente como la conocemos hoy en día. En este ejemplo vamos a deducir de cero la ley de Ohm. Este ejercicio se puede hacer en casa con datos reales si se dispone de un polímetro (dos mejor) y una fuente de alimentación con tensión regulable. Este ejercicio pueden hacerlo niños sin problema.

Olvida todo lo que sepas de la ley de Ohm

Es posible apreciar que en un circuito con una bombilla, si introducimos una pieza cerámica, la intensidad de la bombilla disminuye.

Cuando la corriente no atraviesa la resistencia
Cuando la corriente atraviesa la resistencia

¿Qué ocurre exactamente? ¿Por qué la bombilla tiene menos intensidad en el segundo caso?

Vamos a aislar la resistencia. Ponemos un voltímetro y un amperímetro y vamos cambiando la tensión de entrada. Anotamos la corriente medida en cada caso.

Podemos intentar hacer un ajuste lineal a estos datos. De modo, que una vez sepamos la intensidad, podamos predecir el voltaje.

Como Pearson nos da un número muy próximo a 1, podemos definir un modelo matemático siguiendo la regresión lineal.

Este modelo matemático se define así:

Y es lo que se conoce como la Ley de Ohm. En realidad, la ordenada en el origen tiene un valor muy cercano a cero, podemos quitarlo.

Así nos queda un modelo capaz de predecir lel voltaje en base a la intensidad y la pendiente de la recta. Ahora puedes probar cambiando la resistencia y observando que siempre ocurre lo mismo. Tenemos el voltaje por la pendiente del modelo. Este valor de la pendiente es lo que en física se llama resistencia y se mide en ohmios. Así se nos queda entonces la ley de Ohm que todos conocemos:

\(
V= IR
\)

En este caso la pendiente equivalía a 4.94 Ω, valor muy cercano a los 5 Ω que dice el fabricante.

Autómatas celulares unidimensionales en Python

Estaba yo leyendo este verano un libro titulado Think Complexity cuando en un capítulo empezó a hablar de los autómatas celulares unidimensionales. El tema me interesó y por eso esta entrada. Veamos primero a qué nos referimos cuando hablamos de esto.

Cuando hablamos a autómatas celulares, nos referimos a pequeñas entidades independientes pero que interaccionan entre sí. Celulares porque son la unidad elemental del universo donde van a existir y autómatas porque deciden por ellas mismas, basadas en un conjunto de reglas predefinido, cuando el tiempo avanza de forma discreta (es decir, a pasos).

Este concepto abstracto puede visualizarse con facilidad si nos imaginamos una rejilla. Cada celda es una célula capaz de cambiar su estado según su entorno.

Los autómatas celulares fueron objeto de estudio de Stephen Wolfram, matemático conocido por haber diseñado el programa Mathemathica y Wolfram Alpha.

Los autómatas celulares unidimensionales son aquellos que forman parte de un universo unidimensional. Es decir, cada célula tiene una vecina a su izquierda y a su derecha. En los bordes se pueden definir varios comportamientos pero el concepto no varía. Pensemos en ello como una tira de celdas.

El estudio de estos autómatas es interesante, ya que pueden generarse patrones/situaciones muy complejas en base a unas reglas sencillas.

¿Cómo se definen las reglas?

Wolfram usó un sistema para definir las reglas de estos autómatas que hoy conocemos como Wolfram Code. Se basa en definir una tabla con los estados presentes de la célula y sus vecinas así como el valor que deberá adoptar en esa situación la célula. Como Wolfram usó células con solo dos estados, todo está en binario, y la parte baja de la tabla es un número de 8 bits. Este número se suele pasar a decimal y así se identifica.

Estados presentes 111 110 101 100 011 010 001 000
Estado futuro 0 0 1 1 0 0 1 0

Esta tabla representa la Regla 50, porque 00110010 en binario es 50.

¿Cómo se representan?

Una manera muy interesante de representar estos autómatas es poner cada paso en una fila distinta dentro de una imagen.

Una vez que sabemos esto vamos a hacer un programa en Python que nos permita observar la evolución de estos autómatas.

Usaremos el procedimiento original, que es empezar con todos los estados de los autómatas en 0 salvo el del autómata central, que será 1.

La clase Automata

La clase autómata va a contener las reglas, así como un ID y el estado que posee. Además, por cuestiones técnicas conviene guardar el estado anterior que tuvo.

Como podéis ver, rules no lleva self, es decir, va a ser una variable compartida entre todas las instancias de Automata. Esto es porque las reglas son idénticas a todos los autómatas.

La clase World

Ahora vamos a definir el universo donde residen estos autómatas. Este universo almacena una lista con los autómatas, se encarga de actualizarlos según las normas y de dibujarlos usando PIL. También he insertado el código que codifica las normas según el número en decimal.

Con esto ya lo tenemos casi todo. Ahora faltaría poner en marcha todo. La idea es simplemente crear una instancia de World, hacer unas cuantas llamadas a add, y después ir haciendo el ciclo update/draw_row. Una vez hayamos acabado, hacemos save y obtendremos un PNG con la imagen.

Código completo

Algunas reglas importantes

Regla 30

Una de las más importantes a nivel matemático. Ha sido objeto de mucho estudio, sin embargo no vamos a entrar en detalles más allá de su aspecto visual.

Vista ampliada

Regla 110

Esta regla es también muy interesante. ¡Se demostró que era Turing completa!

Vista en detalle

Regla 126

Esta regla no es tan importante, pero personalmente me parece muy bonita.

Vista ampliada

Reglas 57 y 99

Son dos reglas isomorfas. Es decir, son en realidad la misma regla pero aplicada a lados distintos. Elijo estas dos porque se aprecia muy bien el isomorfismo.

Regla 57
Regla 99

Regla 169

Vista en detalle

Regla 129

Regla 90

Es el famoso triángulo de Sierpinski.

Regla 150

Regla 105

Esta regla no tiene isomorfo.

En este artículo no he querido entrar en las complejidades matemáticas de todo esto. Es algo que todavía no entiendo así que no sería sincero por mi parte exponerlo.

Bonus: Richard Feynman y Steve Jobs

Quien me conoce sabe de sobra que uno de los personajes de la historia que más ha influido en mi vida es Richard Feynman. Debo reconocer que entré en un estado de éxtasis al descubrir que Feynman y Wolfram no solo trabajaron juntos, sino que lo hicieron alrededor de la regla 30 antes mostrada. También me sorprendió que Steve Jobs y Wolfram resultasen ser amigos de toda la vida. No dejo de sorprenderme de los contactos de ciertos personajes históricos entre sí.

Feynman a la izquierda, Wolfram a la derecha

Documentación con rustdoc

La documentación es una parte importante y muchas veces olvidada de un proyecto de software. Rust cuenta desde el principio con una herramienta destinada a que escribir documentación sea poco doloroso, se trata de rustdoc.

Comentarios en Markdown

Los comentarios de rustdoc empiezan con 3 barras y pueden llevar Markdown. Se escriben encima de la función, estructura, etc que describamos.

Mencionar que el código del comentario es sometido a tests también por cargo test. Para generar la documentación basta con escribir cargo doc y Cargo generará la documentación en formato HTML.

Consultando documentación

El mejor sitio para leer la documentación de Rust es Docs.rs. Docs.rs ejecuta cargo doc a todas las crates de Crates.io y las expone al público sin costo.

 

Tests en Rust

Mientras los programas no puedan verificarse de forma matemática de forma sencilla, la única manera de asegurarse que un programa más o menos funciona es con tests. Rust soporta tests de forma nativa. Gracias a la directiva #[test].

Definiendo una función de test

Una función de test es cualquiera que lleve la directiva #[test]. Normalmente, estos tests se dejan dentro de un módulo con la directiva #[cfg(test)].

La macro assert_eq! provocará un panic! si sus argumentos no coinciden. También es posible hacer fallar los tests llamando a panic! manualmente. ¿Cómo se ejecutan estos test te preguntarás? Sencillo, con cargo test. Automáticamente Cargo selecciona las funciones de test y las ejecuta todas, generando un informe de tests existosos y tests que han fallado.

Obviamente, existen más opciones dentro de los tests. assert! comprobará que una expresión sea true. #[should_panic] se deberá indicar en aquellas funciones de test en lo que lo correcto sea que ocurra un panic!. Por otro lado, la trait Debug es interesante.

Es posible ejecutar tests de forma manual, con cargo test NOMBRE_TEST.

 

Cargo y módulos en Rust

El software va creciendo poco a poco en un proyecto y es vital que esté bien organizado. Rust incluye un potente sistema de módulos para manejar estas situaciones y Cargo nos ayuda a gestionar las dependencias del proyecto.

Módulos

Los módulos se definen con la palabra reservada mod y por defecto todo es privado. Es necesario indicar pub para hacer esos campos públicos. Con use podemos acortar la manera con la que hacemos referencia a los componentes del módulo. Así pues, use es similar a using namespace de C++.

Los módulos resultan una opción más interesante si tenemos múltiples archivos. Por norma general, si un módulo se define en el fichero de entrada del compilador, este se busca en dos sitios:

  • Un fichero nombre_del_modulo.rs
  • Un fichero nombre_del_modulo/mod.rs (opción recomendada para módulos de varios archivos)

Por lo que el código anterior puede sustituirse por esto en el fichero principal:

Y esto en un fichero network.rs

Podemos usar use-as para modificar el nombre con el que se importa una función.

 

Crates

El código en Rust se organiza en crates, que son colecciones de módulos que se distribuyen. Piensa en ello como si fuesen gemas de Ruby o librerías de C++. Las crates se pueden compartir con la gente. El mayor repositorio de crates de la comunidad Rust es crates.io.

Las crates se pueden generar con rustc o con cargo, aunque es habitual usar la segunda opción. Vamos a ver primero como usaríamos una crate externa. En este caso voy a usar regex. Sé que en esta parte no vas a saber como descargar regex y usarla, pero voy a explicar el código.

Básicamente se usa extern crate para importar una crate externa, en este caso regex. Con use lo único que hacemos es acortar la manera de llamar a las funciones. Si no estuviese podríamos poner perfectamente regex::Regex::new en vez de Regex::new y funcionaría.

Cargo

Cargo es una herramienta que cumple dos funciones dentro del ecosistema Rust. Por un lado es un gestor de dependencias y por otro lado gestiona también las compilaciones.

En su parte de gestor de dependencias nos permite especificar que crates necesita el proyecto para compilar, su versión y si la crate lo soporta, con qué características activadas.

En su parte de gestor de compilaciones, realiza la compilación con rustc y añade las flags pertinentes al compilador. También se le pueden especificar cosas más complejas, aunque en Rust no suele ser habitual tener configuraciones de compilación excesivamente complejas.

Todo lo relativo a Cargo se define en el archivo Cargo.toml. Que tiene una estructura similar a esta:

En la sección dependencies puedes añadir línea por línea la crate que te haga falta. Consulta la documentación de cada crate para esto, pues puede haber variaciones.

Comandos básicos de Cargo

Crear un proyecto con Cargo

cargo new –bin mi_proyecto # si queremos que sea ejecutable

cargo new mi_crate # si queremos que sea una crate

Compilar y ejecutar

cargo run

cargo build # solo compilar

cargo build –release # compilar en modo Release

Cargo automáticamente se encarga de obtener las dependencias en la fase de compilación.

Instalar aplicaciones

cargo install APLICACION

Muchas herramientas de pueden distribuir a través de Cargo. Por ejemplo, Racer, un programa que sirve de autocompletado para IDEs como Eclipse o Visual Studio.

Ejecutar tests

cargo test

Generar documentación

cargo doc

Plugins de Cargo

Cargo es extensible gracias a plugins. Algunos interesantes son Clippy, cargo-audit, rustfmt,