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.

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.

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.

Estadística en Python: media, mediana, varianza, percentiles (Parte III)

Siguiendo en nuestro camino de la estadística descriptiva, hoy vamos a ver como calcular ciertas medidas relativas a una variable.

Fichero de ejemplo

Medidas de centralización: media, mediana y moda

La media aritmética se define como la suma de N elementos dividida entre N. Se trata una medida bastante conocida entre la gente, aunque tiene el inconveniente de que es muy susceptible a valores extremos.

La mediana es el valor que dentro del conjunto de datos es menor que el 50% de los datos y mayor que el 50% restante.

La moda es el valor más repetido (solo aplicable a variables discretas).

Para calcular estas medidas, simplemente seleccionamos la variable estadística del DataFrame y usamos los métodos mean, median y mode respectivamente.

Ejemplo: Calcula la media, la mediana y la moda de las notas de los alumnos en el examen

Medidas de posición: cuartiles y percentiles

El concepto es igual al de mediana, salvo que aquí la división ya no es en el 50%. El 25% de las observaciones es menor que el primer cuartil. Los cuartiles abarcan el 25%, 50% y 75% de las observaciones. Los percentiles son una generalización con cualquier porcentaje.

Ejemplo: ¿Cuál es la nota que tiene como mínimo el 10% más nota de la clase?

Este enunciado nos pide calcular el percentil 90.

Usamos quantile y el porcentaje.

El resultado es que el 10% con más nota de la clase ha sacado un 8,8 como mínimo. Mencionar que existen distintos tipos de interpolación para este cálculo. En la referencia podemos consultar cual nos conviene más.

Medidas de dispersión: desviación típica, rango, IQR, coeficiente de variación

La desviación típica mide la dispersión de los datos respecto a la media. Se trata de la raíz cuadrada de la varianza, que en sí misma no es una medida de dispersión. Para calcular la desviación típica usamos std y var para la varianza. (ddof=0 es necesario si quieres seguir la definición de desviación típica y varianza de algunas bibliografías, la razón es que hay un parámetro de ajuste que Pandas pone a 1, pero otras librerías ponen a 0). En Excel es la diferencia que hay entre DESVEST.M (ddof=1) y DESVEST.P (ddof=0).

El rango es la diferencia entre el máximo y el mínimo y el rango intercuartílico o IQR es la diferencia entre el tercer y el primer cuartil.

El coeficiente de variación es una medida que sirve para comparar entre dos muestras, cuál varía más y cuál es más estable. Es una simple división, de la desviación típica sobre la media, sin embargo, SciPy nos ofrece una función ya preparada.

Medidas de asimetría

Para saber si los datos estan repartidos de forma simétrica existen varios coeficientes: Pearson, Fisher, Bowley-Yule, etc

Para no liarnos demasiado, podemos usar la función skew de SciPy.

Para valores cercanos a 0, la variable es simétrica. Si es positiva tiene cola a la derecha y si es negativa tiene cola a la izquierda.

Y con esto hemos visto los datos que se pueden extraer de una sola variable.

Estadística en Python: Pandas, NumPy, SciPy (Parte I)

Recientemente SciPy anunció que lanzaba la versión 1.0 tras 16 años en desarrollo. Con motivo de este acontecimiento voy a realizar una serie de tutoriales sobre estadística en Python. Quiero hacer especial énfasis en la resolución de problemas estadísticos, con problemas reales que vamos a ir resolviendo. Usaremos Python 3, NumPy, Pandas, SciPy y Matplotlib entre otras.

¿Por qué Python para estadística?

Python se ha convertido en uno de los lenguajes más usados en la comunidad de data science. Se trata de un lenguaje cómodo para los matemáticos, sobre el que es fácil iterar y cuenta con unas librerías muy maduras.

Instalando el stack estadístico en Python

Yo voy a usar Python 3. Tanto si estamos en Windows, macOS o Linux podemos instalar NumPy, SciPy, Pandas y Matplotlib con este simple comando

¿Para qué sirve cada cosa?

NumPy sirve para realizar cálculos numéricos con matrices de forma sencilla y eficiente y tanto SciPy como Pandas la usan de forma interna. NumPy es la base del stack científico de Python.

SciPy es una colección de módulos dedicados a diversas áreas científicas. Nosotros usaremos principalmente el módulo stats.

Pandas es una librería que permite manipular grandes conjuntos de datos en tablas con facilidad. Además permite importar y exportar esos datos.

Matplotlib permite realizar gráficos y diagramas con facilidad, mostrarlos en pantalla o guardarlos a archivos.

Estadística descriptiva

Vamos a comenzar con la parte de la estadística que trata de darnos un resumen del conjunto de datos.

Para ello vamos a definir el concepto de variable estadística como la magnitud o cualidad de los individuos que queremos medir (estatura, calificaciones en el examen, dinero en la cuenta,…). Las variables pueden ser cualitativas o cuantitativas y dentro de las cuantitativas pueden ser continuas y discretas.

Fichero de ejemplo

En este post voy a usar este fichero para los ejemplos (guardalo como notas.csv).

Cargando datos

Para cargar datos y manipular los datos usamos Pandas. Pandas permite cargar datos de distintos formatos: CSV, JSON, HTML, HDF5, SQL, msgpack, Excel, …

Para estos ejemplos usaremos CSV, valores separados por comas:

df es un objeto de tipo DataFrame. Es la base de Pandas y como veremos, se trata de un tipo de dato muy flexible y muy fácil de usar. Si ahora hacemos print a df obtenemos algo así:

Tabla de frecuencias

Si tenemos una variable discreta o cualitativa una cosa habitual que se suele hacer es construir la tabla de frecuencias. Con ella sabemos cuantas veces se repite el valor de la variable. Para crear tablas de frecuencia usamos crosstab. En index indicamos la variable que queremos contar y en columns especificaos el nombre de la columna de salida. crosstab devuelve otro DataFrame independiente.

Ejemplo: ¿Cuántos alumnos han sacado un 5 en el examen?

Ejemplo: ¿Cuántos alumnos han aprobado (sacar 5 o más)?

En estos ejemplo usamos loc para devolver las filas que cumplan la condición descrita entre corchetes. En el último ejemplo como el resultado son varias filas, nos quedamos con la parte de las frecuencias y sumamos todo, para así obtener el resultado final.

Diagrama de sectores

Una forma sencilla de visualizar datos que ya han sido pasados por la tabla de frecuencias es el diagrama de sectores. Usando Matplotlib podemos generar gráficos en los que podemos personalizar todo, pero cuyo uso básico es extremadamente simple.

Gráfica generada por Matplotlib

Ejemplo: Haz un diagrama de sectores donde se vea claramente el porcentaje de aprobados frente al de suspensos

Destacar que aquí hemos usado np.array para crear un array de NumPy en vez de una lista nativa de Python.

Diagrama de barras

De forma similar es posible generar un diagrama de barras:

Ejemplo: Genere un diagrama de barras de las notas obtenidas por los alumnos en el examen