Estadística en Python: ajustar datos a una distribución (parte VII)

Ya hemos visto con anterioridad que existen modelos que nos hacen la vida más sencilla. Sin embargo, en esos modelos conocíamos ciertos datos de antemano. ¿Qué pasa si tenemos datos y queremos ver si podemos estar ante un modelo de los ya definidos?

Este tipo de ajuste es muy interesante, ya que nos permite saber si los datos en bruto pueden parecerse a los modelos de Normal u otros y aprovecharlo.

Ajuste de datos a una distribución

Para ajustar datos a una distribución, todas las distribuciones continuas de SciPy cuentan con la función fit. Fit nos devuelve los parámetros con los que ajusta nuestros datos al modelo. ¡Ojo, no nos avisa si es un buen ajuste o no!

Ejemplo: Queremos saber si la altura de los hombres adultos del pueblo de Garray sigue una distribución normal. Para ello tomamos una muestra de 80 alturas de hombres adultos en Garray.
Los datos los tenemos en este CSV, cada altura en una línea:

Vamos a ajustarlo.

En este caso nos informa de que estos datos parecen encajar con los de una distribución normal de media 160,37 y desviación típica 17,41.

¿Cómo de bueno es el ajuste? Kolmogorov

Hemos hecho un ajuste, pero no sabemos qué tan bueno es. Existen varios métodos para poner a prueba los ajustes. Existen varios métodos, siendo los más populares Chi-Cuadrado y Kolmogorov-Smirnov.

Chi-Cuadrado no se puede aplicar directamente sobre distribuciones continuas, aún así voy a explicar como se haría. Sin embargo, primero vamos a probar con Kolmogorov-Smirnov, que en SciPy es ktest.

Este test nos devuelve dos valores: D y el p-valor. Yo voy a fijarme en el p-valor. El p-valor es el nivel mínimo de significación para el cual rechazaremos el ajuste. Cuanto más cerca del 1 esté, más confianza hay en el ajuste, cuanto más cerca de cero esté, menos confianza hay en el ajuste. ¿Pero cómo se interpreta exactamente?

Una forma sencillo de entenderlo es hacer obtener el nivel de significación, que es 1 – NivelConfianza/100. Así para una confianza del 95%, este será de 0.05, para una confianza del 99%, el valor se reduce a 0.01. Si el p-valor es inferior a este nivel de significación, el ajuste probablemente no es correcto.

Aquí algún estadístico más serio se tiraría de los pelos, porque en realidad lo único que se puede demostrar con seguridad es en el caso de que no se ajuste. Si salta que se puede ajustar a una normal, simplemente querría decir que no se puede rechazar la hipótesis de partida (que sea normal), pero no confirma que sea normal.

Esto es algo común en el contraste de hipótesis y con los p-valores en particular y ha suscitado crítica, hasta tal punto que en 2016, la American Statistical Association publicó una serie de consejos para tratar con p-valores. El p-valor en este caso solo puede demostrar que no se ajusta a una normal. No obstante, para ir poco a poco, podemos simplificar esto un poco.

Chi-Cuadrado

Para hacer el test de Chi-Cuadrado primero hay que dividir los datos continuos. Para uso podemos usar la función de NumPy, histogram o la de SciPy binned_statistic. Tenemos que construir una lista con las frecuencias esperadas si siguiéramos al 100% la distribución a la que queremos ajustarnos.

Una vez tengamos una lista con los valores esperados, simplemente podemos comparar las frecuencias reales con las esperadas con chisquare. Esto nos devolverá C y el p-valor. El significado del p-valor es exactamente igual.

Con esto ya hemos visto como empezar con estadística en Python. No sé si haré más artículos de este estilo, si no es así, espero que os hayan gustado esta serie de artículos. El mundo del data science es inmenso y os invito a seguirlo explorando.

Yew, crea webapps al estilo Angular/React en Rust

Hoy os traigo una librería muy potente y muy útil, ahora que Rust compila a WebAssembly de forma nativa. Se trata de Yew, una librería para diseñar frontends, single-page-applications y en general, una librería que es una alternativa a Angular, React, Vue y Elm.

En particular, Yew se basa en The Elm Architecture, por lo que los usuarios de Elm serán los que encuentren más familiar a Yew.

Yew significa tejo. He aquí uno de los tejos más famosos de Valladolid, en la plaza del Viejo Coso

Instalar cargo-web y el soporte a WebAssembly

Antes de empezar a programar necesitaremos instalar un plugin para Cargo, llamado cargo-web, que nos ayudará en el desarrollo web con Rust. Por otro lado, hace falta instalar el soporte de Rust a WebAssembly. Existen tres opciones actualmente: asmjs-unknown-emscripten, wasm32-unknown-emscripten y wasm32-unknown-unknown. Para los primeras opciones hace falta tener instalado Emscripten. Para la última, no hace falta nada, por lo que es mi opción recomendada. Por otro lado, wasm32 significa WebAssembly y asmjs es asm.js, que es simplemente JavaScript y será compatible con más navegadores.

The Elm Architecture

Los principios de The Elm Architecture se basan en: modelo, mensajes, actualización y vista.

Para este tutorial vamos a hacer la típica aplicación de una lista donde guardamos notas. Nuestro modelo se va a componer de una lista de tareas, para simplificar, pongamos que una tarea es simplemente un String y un ID. Entonces también nos hará falta almacenar un contador para ir poniendo IDs. También, al tener un campo de texto, nos hará falta una variable para ir almacenando temporalmente lo que escribe el usuario.

Lo siguiente es diseñar los mensajes. Los mensajes interactúan con el modelo y desencadenan una actualización de la vista. En esta aplicación solo nos hacen falta dos mensajes: añadir mensaje y borrar mensaje. Pero como tenemos un campo de texto, tenemos que introducir un mensaje Change y siempre viene bien un mensaje que no haga nada.

Una vez hecho esto pasamos a crear la función update, en la que hacemos distintas cosas según el mensaje que recibamos.

Así pues, si se lanza un mensaje Msg::Add lo que hacemos es copiar el valor de la variable temporal input, crear una nueva tarea con su ID y añadirla a la lista de tareas. Ya está. Yew mágicamente actualizará la página para reflejar que la lista de tareas ha sido modificada. Lo mismo pasa con Remove.

Ahora vamos a las vistas. Una vista es una función que devuelve Html<Msg> y se pueden componer varias funciones así. En nuestro caso, tenemos una vista principal donde se ve un campo de texto y un sitio donde se ejecuta un bucle for con las tareas del modelo. Y a cada tarea se le aplica la vista view_task.

La macro html! nos permite escribir HTML directamente en Rust, con algunas diferencias (¡prestad atención a las comas!). También nos permite introducir código Rust (entre llaves) y closures (observad onclick, oninput y onkeypress).

Finalmente en el método main, inicializamos el modelo y llamamos a program, que empieza a ejecutar Yew.

El código final queda así.

Ejecutando la aplicación web

Usando cargo web, es muy sencillo generar la aplicación web. Simplemente ejecuta:

El resultado, se montará en un mini servidor web. Si accedes a la URL que indica Cargo con tu navegador web, verás algo similar a esto:

Añade items y bórralos. Observa como la aplicación funciona perfectamente.

Distribuyendo la aplicación

Realmente cargo web ha hecho muchas cosas por nosotros. Si nosotros queremos usar Yew en la vida real, no usaremos cargo web. Para ello, compilamos la aplicación web:

Y accedemos a la carpeta target/wasm32-unknown-unknown/release. Allí encontraremos dos archivos que vamos a necesitar. Uno acabado en .js y otro acabado en .wasm. Ambos ficheros deberemos copiarlos donde queramos usarlos. Por último, será necesario un fichero HTML. En el fichero HTML solo hace falta cargar el fichero JS. Yew hará el resto.

Si quieres saber más, puedes descargarte el código de ejemplo que se encuentra en GitHub.

Leer y escribir JSON en Rust con Serde

JSON es un formato muy popular para guardar y transmitir información. En Rust podemos leer y escribir JSON de forma transparente gracias a Serde. Serde es una librería para Rust, que permite transformar las structs nativas de Rust en ficheros JSON, YAML, BSON, MsgPack, XML y viceversa. Serde está compuesto de varios plugins, uno para cada formato de archivo, así que vamos a necesitar el plugin de JSON. Es lo que se dice serializar (guardar) y deserializar (leer).

Deserializando JSON con Serde

Para deserializar necesitamos primero definir el struct en Rust que se corresponde al JSON que vamos a cargar. Además, hay que incluir una cláusula derive e indicar que derivamos Deserializable.

Una vez hecho eso podemos realizar la transformación con Serde. Serde admite varios métodos de entrada. En mi caso voy a usar from_reader, que permite transformar desde cualquier std::io::Read. También se podría usar from_str, que permite leer desde un String.

Así pues un programa que lea Landmarks de un fichero JSON como este:

Quedaría así:

Nótese como usamos Vec<Landmark> para indicar que vamos a cargar un JSON que es un array de Landmarks.

Serializar JSON

Ahora si queremos generar un archivo JSON es muy sencillo. Simplemente marcamos Serialize en la misma estructura que queramos serializar. Posteriormente en Serde podemos utilizar diferentes métodos, como to_writer o to_string.

Personalizando la serialización y la deserialización

Serde permite a través de atributos, definir de forma precisa que ocurre con los elementos. En caso de no poner nada, por ejemplo, Serde usará los nombres del struct de Rust en los ficheros y si hay elementos de más, en el fichero al leer, los ignorará.

Por ejemplo, si queremos que en el fichero JSON, en vez de tener name sea nombre, podemos usar rename.

Podemos también no querer guardar algún elemento del struct, usamos skip.


Si queremos que Serde falle en caso de que en el JSON haya más campos de los que hemos definido, usamos deny_unkown_fields.

Y si queremos que el nombre del struct sea distinto al que crea oportuno Serde, podemos redefinirlo también:

Por último, mencionar que Serde permite serializar/deserializar cosas más complejas, con otros structs, con vectores, con HashMap, entre sus elementos,… Lo único que tendrá que pasar es que esas estructuras a su vez sean serializables/deserializables con Serde (es decir, hayan puesto derive(Serialize,Deserialize)).

¡Feliz navidad y próspero 2018! (con fractales)

Feliz Navidad! Si estás leyendo esto, darte las gracias por seguir leyendo este blog. Me alegro mucho de tener cada vez más lectores y es algo que me anima a escribir más y más. En este especial navideño disfrutaremos de éxitos musicales de la navidad, veremos como hacer un fractal muy navideño en Rust y os pediré vuestra opinión.

Comenzamos con In Dulci Jubilo y este cover de Giuliano Ferrace Leanza

El fractal de Koch

Helge von Koch fue un matemático sueco que en 1904 publica un artículo sobre la que desde entonces se llamaría curva de Koch. Se trata de una curva, generada con unas reglas muy simples. La idea fundamental es que se parte de una recta, que se divide en tres trozos iguales. El trozo del medio se sustituye por un triángulo equilátero. Entonces la curva recorre el primer trozo, sube el triángulo, baja el triángulo y continúa recto.

Si ahora en cada uno de los cuatro trozos rectos aplicamos la misma curva de Koch obtenemos esto:

Según la terminología de Mandelbrot, este fractal es un teragón. Si ahora en vez de aplicar esto sobre una recta, lo hacemos sobre un triángulo equilátero con sus tres trozos rectos iguales, ¿qué obtenemos?

Turtle, una librería de Rust para gráficos tortuga

Quizá entre mis lectores haya alguno que aprendió a programar con LOGO. Una de las características de LOGO era que disponía de una tortuga con la que íbamos dibujando según nos movíamos. Esta manera de dibujar, no es óptima en rendimiento puro, pero es ideal para fractales y para enseñar a programar. Python incorpora en su librería estándar el módulo turtle y en Rust está a punto de salir una librería que soporta una API similar llamada turtle.

Lo primero que tenemos que hacer es definir la curva de Koch con turtle.

Con esto nos valdría, pero no es recursivo. Necesitamos poder aplicar la curva de Koch en nuestras rectas (donde hacemos forward) de forma recursiva. Pero si lo hacemos recursivo de forma infinita no acabará nunca, es por ello que tenemos que indicar cuantos niveles de recursividad queremos.

Esto ya tiene más sentido. Ahora juntemos todo lo necesario:

Si compilamos y ejecutamos esto con Cargo:

Y aquí tenemos al copo de nieve de Koch.

Tenéis el proyecto completo en GitHub: https://github.com/aarroyoc/fractal_koch_rust

Fractal (bis)

Otro fractal muy interesante que programé mientras estaba haciendo el copo de nieve fue este otro. Desconozco si tiene nombre. Es algo más complejo de implementar y lo de los colores me dio bastantes dolores de cabeza hasta que encontré una progresión bella.

Su código es (en este caso Python):

 

Vuestra opinión

Ahora os pido que me déis vuestra opinión sobre el blog. Usad los comentarios de debajo y contadme: ¿Cuál es el mayor problema del blog?, ¿crees que los contenidos salen con frecuencia suficiente y necesaria?, ¿los temas son interesantes?, ¿me voy mucho por las ramas?, … Todo esto lo tendré en cuenta de cara al año que viene.

Aprovecho para recordar que podéis suscribiros a la lista de correo para que os llegue un nuevo correo con cada post, podeís suscribiros al blog por RSS, hay un canal de Telegram y una página en Facebook. En Twitter, GNU Social e Instagram solo tengo perfiles personales, pero si os interesa, allí estoy. Por último en Google+ todavía sigo mandando los artículos.

Si os parece que me merezco una caña, siempre puedes donar usando PayPal, criptodivisas, tu apartamento en Comillas, …

Sin más, feliz próspero año 2018

 

Estadística en Python: distribución binomial, normal y de Poisson (Parte VI)

Después del rollo del capítulo anterior, vamos a entrar en algo más práctico, modelos de probabilidad. Hemos hablado de las funciones de probabilidad y de densidad, sin embargo, nos falta algo muy importante. ¿Cuáles son esas funciones exactamente? Es donde entran los modelos de probabilidad, modelos estadísticos que se pueden ajustar a una variable aleatoria con mejor o peor precisión y que nos dan los valores de la probabilidad. Empecemos, pero antes hagamos un apunte sobre las equivalencias en SciPy:

  • cdf(x) – Función de distribución F(X)
  • sf(x) = 1 – cdf(x)
  • pmf(x) – Función de probabilidad f(x) (distribuciones discretas)
  • pdf(x) – Función de densidad f(x) (distribuciones continuas)
  • ppf(x) – Función inversa a cdf(x). Nos permite obtener el valor correspondiente a una probabilidad.

Distribución Binomial

Un ensayo de Bernouilli se define como un experimento donde puede darse un éxito o fracaso y donde cada ensayo es independiente del anterior. Por ejemplo, un ensayo de Bernoulli de parámetro 0.5 sería lanzar una moneda a cara o cruz (mitad de posibilidades de cara, mitad de posibilidades de cruz).

Si repetimos N veces los ensayos de Bernouilli tenemos una distribución binomial.

\(
X \rightarrow B(N,P)
\)

SciPy nos permite usar binom para trabajar con distribuciones binomiales.
Ejemplo: Un proveedor de DVDs regrabables afirma que solamente el 4 % de los
artículos suministrados son defectuosos. Si un cliente compra un lote de 25
DVDs, ¿cuál es el número esperado de DVDs defectuosos en el lote? Si el cliente
encuentra que 4 de los DVDs comprados son defectuosos, ¿debe dudar de la
afirmación del vendedor?
El número de DVDs defectuosos esperados es el equivalente a decir el número medio de DVDs defectuosos.

Es decir, de media habría 1 DVD defectuoso en el paquete. mean calcula la media de la distribución.

Para saber si hay que fiarse del vendedor vamos a calcular cuál era la probabilidad de que nos tocasen 4 DVDs defectuosos.

Es decir, la probabilidad que ocurriese era del 1%. Podemos sospechar del fabricante. cdf calcula las probabilidades acumuladas. En este caso tenemos que calcular la probabilidad de que hubiese 4 o más fallos, Pr{X>=4}. Una manera fácil de calcularlo es hacer 1-Pr{X<4}. cdf(n) nos permite calcular probabilidades acumuladas hasta N. Otra opción sería simplemente obtener la probabilidad de 0 DVDs defectuosos, 1 DVD defectuoso, de 2 DVDs defectuosos, de 3 DVDs defectuosos, sumarlo y restarlo de 1.

pmf(n) devuelve la probabilidad de que X=N, Pr{X=N} Esto solo tiene sentido en ciertas distribuciones, las discretas, como es el caso de la binomial.

Podemos calcular la gráfica de esta distribución binomial.

En el gráfico también se puede ver que las probabilidades de tener 4 o más DVDs defectuosos son mínimas.

Distribución hipergeométrica

La distribución hipergeométrica es un modelo en el que se considera una población finita de tamaño N en la cual hay M individuos con una determinada característica y se seleccionan n y queremos saber la probabilidad de que haya cierto número de individuos con esa característica en la selección. Para trabajar con estas distribuciones, SciPy trae hypergeom.

Ejemplo: Se formó un jurado de 6 personas de un grupo de 20 posibles miembros de los cuales 8 eran mujeres y 12 hombres. El jurado seelecionó aleatoriamente, pero solamente tenía 1 mujer. ¿Hay motivos para dudarde la aletoriedad de la selección?

La probabilidad de que ocurriese lo que ocurrió es del 18,7%, una probabilidad suficientemente alta como para pensar que no hubo manipulación. Podemos dibujar esta hipergeométrica:

Como se puede observar en la gráfica el caso más probable, con cerca del 35% de posibilidades era que hubiese dos mujeres seleccionadas. Destacar también, que la probabilidad de que haya 7 mujeres en el jurado es cero, porque solo hay 6 plazas en el jurado.

Distribución de Poisson

La distribución de Poisson recoge sucesos independientes que ocurren en un soporte continuo. El número medio de sucesos por unidad de soporte se le conoce como λ y caracteriza la distribución. poisson nos permite crear distribuciones de este tipo.

Algunos ejemplos de distribuciones de Poisson: número de clientes que llegan cada hora a cierto puesto de servicio, número de averías diarias de un sistema informático, número de vehículos que pasan diariamente por un túnel, número de defectos por kilómetro de cable, …

Ejemplo: La impresora de una pequeña red informática recibe una media de 0.1 peticiones por segundo. Suponiendo que las peticiones a dicha impresora son independientes y a ritmo constante, ¿cuál es la probabilidad de un máximo de 2 peticiones en un segundo? Si la cola de la impresora tiene un comportamiento deficiente cuando recibe más de 10 peticiones en un minuto, ¿cuál es la probabilidad de que ocurra esto?

Variable Y: número de peticiones a la impresora en un minuto (y la probabilidad de que suceda)

Distribución exponencial

Para modelizar el intervalo entre dos sucesos consecutivos que siguen una distribución de Poisson se usa la distribución exponencial de parámetro λ.

Ejemplo: El proceso de accesos a una página web se produce de una forma estable e independiente, siendo el intervalo entre dos accesos consecutivos una v.a. exponencial. Sabiendo que, de media, se produce un acceso cada minuto,¿cuál es la probabilidad de que no se produzcan accesos en 4 minutos? y ¿cuál esla probabilidad de que el tiempo transcurrido entre dos accesos consecutivos sea inferior a 90 segundos?

Esta distribución en SciPy es un poco rara, ya que no está implementada como podría esperarse.

Distribución normal

Probablemente el modelo de distribución más usado y conocido. Lo usamos para describir variables reales continuas.

Ejemplo: La duración de un determinado componente electrónico, en horas, es una v.a. que se distribuye según una N(2000,40). ¿Cuál es la probabilidad de que la duración de una de esas componentes sea superior a 1900 horas? ¿y de que esté entre 1850 y 1950 horas?

Podemos representar esta variable.

Estos modelos no son perfectos, pero son lo suficientemente flexibles para ser un buen punto de partida.