La Paradoja de Simpson

Imagen de John Hain en Pixabay 

Introducción

Por casualidad me encontré con la llamada Paradoja de Simpson, y desde el primer momento, y como buena paradoja, me llamó la atención. Me propuse escribir una entrada en el blog sobre ella. Buscando por la web encontré muchísimos artículos al respecto, de los que destaco varios en Medium, aquí y aquí por ejemplo, todos muy buenos y explicativos. Y en uno de ellos conocí al autor Judea Pearl….y vaya descubrimiento. Lo que, a priori, me pareció un entretenido juego de números explotó y me permitió descubrir la disciplina (¡o ciencia!) de la «causal inference». Empecé leyendo (y tratando de entenderlos en su totalidad!) sus papers, aquí y aquí, el segundo bastante más técnico al basarse en su libro (Causality: Models, Reasoning and Inferencies, 2000) y sobre todo el capítulo 6 de su último libro (The Book of Why, 2019). Libro que comprende una explicación muy asequible de los modelos que el autor ha ido desarrollando durante su carrera. Sin haber leído estas referencias no podríamos más que hablar de un números y las curiosidades de las fracciones resultantes de suma de otras, pero ahora podemos enfocar la explicación, aún siendo simple, de una manera más «científica»…vamos con ello.

Leer más »

Gestión de followers de Twitter con Python

Seguimiento de Followers

Existen varias herramientas para gestionar cuentas de twitter que te permiten controlar el movimiento de seguidores: quien te sigue, te deja de seguir, etc…
Para practicar un poco construyendo apps me he puesto a codificar un sencillo programa que te gestione los followers de tu cuenta sin necesidad de instalar las herramientas mencionadas. Este código puede ser el embrión sobre el que añadir más funcionalidades. Ahí queda trabajo pendiente.

Leer más »

El Sesgo en los algoritmos de decisión automáticos

Resumen

Es inherente a los algoritmos de IA sufrir un sesgo que llevado a los extremos puede reforzar la discriminación en determinados entornos, y forzar decisiones injustas. Explicaremos la base teórica, en un modelo muy simplificado, de este sesgo y añadiremos algunos ejemplos reales. Y veremos como se intentar evitar tanto desde el punto de vista legal como desde el punto de vista más técnico. Tanto los organismos responsables de la UE como muchas de las empresas del sector están lanzando recomendaciones para el diseño y operación de los algoritmos de decisión. De estas expondremos las más relevantes, que nos han de servir como guía a la hora de construirlos.

(Imagen de Somchai Chitprathak en Pixabay)

Leer más »

Teoremas «límite» de la Estadística: Ley de los Grandes Números (LLN: Large Numbers Law) (1)

Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
view raw lln.ipynb hosted with ❤ by GitHub

Pandas

Introducción al Dataframe de Pandas con Python (I)

En esta entrada vamos a trabajar con Pandas, una estructura de datos fundamental para trabajar organizar los datasets en DataScience. Es una libreria muy completa, muy documentada en multitud de blogs y que veremos como nos facilita el trabajo. Empecemos!:

import pandas as pd
import numpy as np
Leer más »

Gestión de ‘listas’ en Python

Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
view raw list.ipynb hosted with ❤ by GitHub

Análisis de tweets en python (II)

  

Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

  

Análisis de tweets con Python (I)

Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

Análisis del fichero de Licencias de Taxis de Madrid

Análisis de la planta de Taxis de Madrid

Vamos a analizar el fichero de «Taxis»del Ayuntamiento de Madrid, con información sacada del portal de OpenData : https://datos.madrid.es/portal/site/egob/menuitem.c05c1f754a33a9fbe4b2e4b284f1a5a0/?vgnextoid=4f16216612d39410VgnVCM2000000c205a0aRCRD&vgnextchannel=374512b9ace9f310VgnVCM100000171f5a0aRCRD

Como siempre importamos las librerias necesarias : pandas, numpy, matplotlib,datetime..y añadimos calendar. Usaremos esta última para sacar el literal de un mes según su nº de orden.

import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import datetime
import matplotlib.dates as mdates
%matplotlib inline
import matplotlib.ticker as mtick
from matplotlib.ticker import FuncFormatter
from mpl_toolkits.mplot3d import Axes3D
import calendar
pd.options.display.float_format = '{:,.1f}'.format

Preparamos una texto para incluirlo en cada gráfico como fuente…

fuente='Fuente : Ayuntamiento de Madrid, http://datos.madrid.es'

Construimos la URL de la fuente de datos :

path_web='https://datos.madrid.es/egob/catalogo/300219-0-taxi-flota.csv'

Leemos los datos desde su localizacion en ‘path_web’, en este fichero tenemos los datos de Junio de 2018.

taxis=pd.read_csv(path_web,sep=";",encoding='windows-1250',index_col=False)

Vemos unos ejemplos de lineas para asegurar que todo ha ido bien. Podemos chequear también las columnas

len(taxis)

15652

taxis.describe()
CĂłdigo Cilindrada
count 15,652.0 15,652.0
mean 1,348,913.2 1,687.2
std 100,096.0 171.3
min 550,016.0 0.0
25% 1,269,989.8 1,598.0
50% 1,346,296.5 1,598.0
75% 1,433,151.5 1,798.0
max 1,514,196.0 3,498.0

«`python
taxis.columns
«`

Index([‘Código’, ‘Matrícula’, ‘Fecha Matriculación’, ‘Marca’, ‘Modelo’, ‘Tipo’,
‘Variante’, ‘Clasificación medioambiental’, ‘Combustible’, ‘Cilindrada’,
‘Potencia’, ‘Número de Plazas’,
‘Fecha inicio de prestación del servicio de taxi’, ‘Eurotaxi’,
‘Régimen Especial de Eurotaxi’,
‘Fecha inicio Régimen Especial Eurotaxi’,
‘Fecha fin Régimen Especial Eurotaxi’, ‘Fecha’],
dtype=’object’)

Nos aseguraremos de que la columna ‘Fecha Matriculacíon’ tiene el formato Datetime para trabajos posteriores, y la ordenamos a continuación

taxis['Fecha Matriculación']=pd.to_datetime(taxis['Fecha Matriculación'],format='%d/%m/%Y')
taxis=taxis.sort_values(['Fecha Matriculación'])

Busquemos el Taxi con la matrícula más antigua :

print ('Matricula {} de fecha {}'.format(taxis['Matrícula'].iloc[0],taxis['Fecha Matriculación'].iloc[0].strftime('%d/%m/%Y')))

Matricula 2239DTJ de fecha 13/12/2005

..y con la matrícula más moderna :

print ('Matricula {} de fecha {}'.format(taxis['Matrícula'].iloc[-1],taxis['Fecha Matriculación'].iloc[-1].strftime('%d/%m/%Y')))

Matricula 8333KLS de fecha 12/06/2018

Calculemos la lista de número de Taxis por Marca

fig = plt.figure(1, (6,4))
ax = fig.add_subplot(1,1,1)
ax=taxis['Marca'].value_counts().plot.bar()
ax.locator_params(axis='y',nbins=10)
ax.set_xlabel('Marca')
ax.set_ylabel('Número de coches',size=16)
ax.grid(axis='y')
ax.set_title('Nº de Taxis por Marca')
fig.suptitle(fuente,size=10,x=0,y=-.4)
fig.savefig('Taxis por Marca',bbox_inches = 'tight')

output_32_0

..ahora ordenaremos por combustible usado

fig = plt.figure(1, (6,4))
ax = fig.add_subplot(1,1,1)
ax=taxis['Combustible'].value_counts().plot.bar()
ax.locator_params(axis='y',nbins=10)
ax.set_xlabel('Combustible')
ax.set_ylabel('Número de coches',size=16)
ax.grid(axis='y')
ax.set_title('Nº de Taxis por Combustible')
fig.suptitle(fuente,size=10,x=0,y=-.4)
fig.savefig('Taxis por Marca',bbox_inches = 'tight')

output_34_0

Ahora trabajemos estos dos últimos campos de manera conjunta. Empezamos con un ‘groupby’ por estos dos campos, añadiendo ‘size()’ para que nos calcule el tamaño de cada campo :

tabla=taxis.groupby([taxis['Marca'],taxis['Combustible']]).size()
tabla

Marca Combustible
CHEVROLET DIESEL 7
GASOLINA TRANSFORMADO GLP 1
CITROEN DIESEL 340
GLP / GASOLINA 345
DACIA DIESEL 281
GLP / GASOLINA 896
FIAT DIESEL 54
GLP / GASOLINA 276
FORD DIESEL 145
GASOLINA-ELECTRICIDAD 4
HYUNDAI DIESEL 7
MERCEDES-BENZ DIESEL 199
GASOLINA – GAS NATURAL 9
GASOLINA TRANSFORMADO GLP 1
NISSAN ELECTRICO 13
OPEL DIESEL 1
PEUGEOT DIESEL 992
PEUGEOT DIESEL 1
RENAULT DIESEL 134
SEAT DIESEL 1875
GASOLINA – GAS NATURAL 199
GASOLINA TRANSFORMADO GLP 48
GLP / GASOLINA 612
SKODA DIESEL 4412
GASOLINA TRANSFORMADO GLP 14
GLP / GASOLINA 207
SSANGYONG GASOLINA 6
GASOLINA TRANSFORMADO GLP 1
GLP / GASOLINA 6
TOYOTA GASOLINA-ELECTRICIDAD 4067
VOLKSWAGEN DIESEL 478
GASOLINA – GAS NATURAL 21
dtype: int64

Ya tenemos la tabla lista, ahora vamos a representarla. He encontrado este modelo de mapa de temperatura de SeaBorn que queda muy bien :

tabla=tabla.reset_index()
import seaborn as sns
sns.set()
matriz=tabla.pivot('Marca','Combustible',0)
f, ax = plt.subplots(figsize=(15, 6))
sns.heatmap(matriz, fmt=".0f", annot=True,linewidths=1, ax=ax, cmap="YlGnBu")

output_39_1

Insertemos ahora una columna con la antigüedad de cada Matricula, simplemente calculando la diferencia entre ‘ahora’ y la fecha de matriculación, que vendrá en días (.days) y la pasamos a años equivalentes dividiendo por 365 (nos olvidamos de los bisiestos!)

ahora=datetime.datetime.today()
taxis['antiguedad']=taxis['Fecha Matriculación'].apply(lambda x : (ahora-x).days/365)

A continuación representamos una gráfica con el nº de matriculaciones ordenadas por meses :

fig = plt.figure(1, (15,6))
ax = fig.add_subplot(1,1,1)

data_series=taxis.groupby([taxis['Fecha Matriculación'].dt.year,taxis['Fecha Matriculación'].dt.month]).size()

ax.set_xlabel('Marca')
ax.set_ylabel('Número de coches',size=16)

ax.grid(axis='y')
ax.set_title('Nº de Matriculaciones de Taxis por mes')

bar_width=0.8
opacity=1
error_config = {'ecolor': '0.3'}
rects1 = ax.bar(np.arange(len(data_series.index.tolist())), data_series.values,bar_width,alpha=opacity)
def format_func(value, tick_number):
anno,mes=data_series.index.tolist()[int(value)]
mes=calendar.month_abbr[mes]
return mes+'-'+str(anno)[-2:]

ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))

# Preparamos una lista con los tickts que nos interesan : el primero, el último y todos los meses Enero y Julio
lista_ticks=[]
lista_ticks=[i for i in np.arange(3,len(data_series.index.tolist())) if data_series.index[i][1] in (1,7)]
lista_ticks.extend([0,len(data_series)-1]) #añade el primero y el último
ax.xaxis.set_ticks(lista_ticks)
fig.suptitle(fuente,size=15,x=0.5,y=0)
fig.savefig('matriculas por mes',bbox_inches = 'tight')
plt.show()

output_43_0

Utilicemos el campo ‘antiguedad’ para un par de cálculos y para generar una gráfica :

'Antiguedad media del parque de taxis : {:.1f} años, con desviación +/-{:.1f} años'.format(taxis['antiguedad'].mean(),taxis['antiguedad'].std())

‘Antiguedad media del parque de taxis : 4.0 años, con desviación +/-2.4 años’

Se nos ocurre generar una gráfica con la antigüedad media por Marca :

taxis.groupby(['Marca'])['antiguedad'].describe().sort_values(by='mean',ascending=False)
count mean std min 25% 50% 75% max
Marca
HYUNDAI 7.0 8.7 0.5 8.2 8.4 8.8 8.9 9.5
OPEL 1.0 8.7 nan 8.7 8.7 8.7 8.7 8.7
CHEVROLET 8.0 8.1 0.3 7.7 7.8 8.0 8.4 8.5
PEUGEOT 1.0 6.4 nan 6.4 6.4 6.4 6.4 6.4
SKODA 4,633.0 5.0 2.4 0.5 2.9 5.0 6.8 12.6
RENAULT 134.0 4.9 1.8 0.2 4.4 5.2 6.2 8.8
TOYOTA 4,067.0 4.4 2.5 0.1 2.2 4.9 6.5 10.4
MERCEDES-BENZ 209.0 4.1 2.3 0.2 2.7 3.6 5.1 10.0
SEAT 2,734.0 3.9 2.0 0.1 2.4 4.2 5.1 11.0
VOLKSWAGEN 499.0 3.7 2.6 0.1 0.9 3.6 4.9 10.0
PEUGEOT 992.0 3.3 1.7 0.5 2.2 3.1 4.2 9.5
DACIA 1,177.0 2.4 1.0 0.1 1.5 2.8 3.2 3.8
NISSAN 13.0 2.1 0.7 1.0 2.0 2.0 2.2 3.8
CITROEN 685.0 2.0 1.1 0.1 1.2 2.0 2.7 9.3
FIAT 330.0 0.8 1.2 0.1 0.3 0.5 0.9 9.7
FORD 149.0 0.8 0.5 0.1 0.5 0.6 0.9 2.6
SSANGYONG 13.0 0.4 0.2 0.1 0.2 0.5 0.6 0.7

«`python
taxis.groupby([‘Marca’])[‘antiguedad’].mean().sort_values(ascending=False).plot.bar()
«`

output_48_1

Y jugando un poco con Seaborn hemos encontrado este formato de gráficas que muestra muy claramente la contribución de cada generación de coches a esas medias :

sns.set(style="whitegrid", palette="muted")
n_de_marcas=8
#filtramos a las primeras n_marcas por volumen
#datos=taxis.loc[taxis['Marca'].isin(taxis_marca[:n_de_marcas].index)]
lista_topn_marcas=taxis['Marca'].value_counts().index[:n_de_marcas]
datos=taxis.loc[taxis['Marca'].isin(lista_topn_marcas)]
# Draw a categorical scatterplot to show each observation
f, ax = plt.subplots(figsize=(15, 6))
sns.swarmplot(x="Marca", y="antiguedad", ax=ax, data=datos)

output_50_1

Vemos que se ha dejado de vender Skodas. Hay varios huecos en Volkswagen. Fiat ha empezado a vender en los dos últimos años. Peugeot ha dejado de vender también, y Dacia ha recupeado tras algún mes sin vender.

sns.set(style="whitegrid", palette="muted")
n_de_marcas=2
#filtramos a las primeras n_marcas por volumen
#datos=taxis.loc[taxis['Marca'].isin(taxis_marca[:n_de_marcas].index)]
lista_topn_marcas=taxis['Marca'].value_counts().index[:n_de_marcas]
datos=taxis.loc[taxis['Marca'].isin(lista_topn_marcas)]
# Draw a categorical scatterplot to show each observation

f, ax = plt.subplots(figsize=(15, 6))

sns.swarmplot(x="Marca", y="antiguedad",hue='Modelo', ax=ax, data=datos)
#lg=sns.legend_
plt.legend(ncol=2, loc='upper left')

output_52_1

Eso es todo!

Paradoja del análisis médico generalizado

Problema del Análisis médico

Una determinada prueba clínica para detectar una enfermedad tiene una sensibilidad del 98%, es decir : en el caso de personas afectadas acierta en el 98% de los casos. El el 2% restante da resultado negativo, cuando debería de ser positivo (Falso Negativo).
Y tiene una especificidad del 99.0%, es decir : en el caso de personas sanas acierta en el 99.0% de los casos, y yerra en el 1.0% restante, dando un resultado positivo cuando en realidad el paciente está sano (Falso Positivo). La tasa de enfermos en la población del 0.5%.
La pregunta es : si hacemos un análisis a una persona elegida al azar de entre toda la población, y el resultado da positivo : ¿qué posibilidades hay de que realmente esté enferma?

Solución

La teoria de Bayes dice :
p(B_{i}| A)=\frac{p(A | B_{i}) * p(B_{i})}{\sum_{k=1}^{n}p(A | B_{k})*p(B_{k})}

Vamos a aplicarlo con la siguientes nomenclatura:
A_{+} : posibilidad de que el análisis de positivo
A_{-} : posibilidad de que el análisis de negativo
E_{+} : la persona está enferma
E_{-} : la persona está sana

Aplicando los datos que nos han pasado en el enunciado :

p(A_{+}| E_{+})=0.98\\ p(A_{-}| E_{+})=0.02\\ p(A_{+}| E_{-})=0.01\\ p(A_{-}| E_{-})=0.99

p(E_{+})=0.005

El dato que nos piden es : p(E_{+} | A_{+}) probabilidad de que esté enfermo (E_{+}) habiendo dado positivo en el análisis (A_{+})

p(E_{+}| A_{+})=\frac{p(A_{+} | E_{+}) * p(E_{+})}{\sum_{k=1}^{2}p(A_{+} |t E_{k})*p(E_{k})}

desarrollando el $\sum$ :
p(E_{+}| A_{+})=\frac{p(A_{+} | E_{+}) * p(E_{+})} {p(A_{+} | E_{+})*p(E_{+})+p(A_{+} | E_{-})*p(E_{-})}

Usemos Python para el cálculo:

Empecemos importando una serie de librerías estandares en estas entradas :

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import re
import math

Criterio : mayúscula es positiva, minúscula es negativa

pAE=0.98
paE=1-pAE
pAe=0.01
pae=1-pAe
pE=0.005
pe=1-pE
pEA=pAE*pE/(pAE*pE+pAe*pe)
'La probabilidad de que esté enfermo dando positivo en el test es del {:1.2f}%'.format(100*pEA)
'La probabilidad de que esté enfermo dando positivo en el test es del 33.00%'

Esa es la razón por la que no tienen sentido análisis a la población en general cuando la tasa de afectados es muy baja.

Para entender el resultado hagamos una extrapolación para España (nº aprox)
Población : 46,5 millones, que podemos dividir en :

poblacion=46500000
poblacion_e=poblacion*pe
poblacion_E=poblacion-poblacion_e
'Poblacion afectada {:,.0f}'.format(poblacion_E)
'Poblacion afectada 232,500'
'Poblacion NO afectada {:,.0f}'.format(poblacion_e)
'Poblacion NO afectada 46,267,500'

aplicando p(A|e) a la poblacion no enferma nos encontrariamos con :

poblacionAe=poblacion_e*pAe
poblacionAe
462675.0
'es decir : tendriamos {:,.0f} False Positivos'.format(poblacionAe)
'es decir : tendriamos 462,675 False Positivos'

Y si hacemos el test solo a la poblacion afectada : cuantos positivos (verdaderos en este caso) tendriamos?

poblacionAE=poblacion_E*pAE
'{:,.0f} True Positivos'.format(poblacionAE)
'227,850 True Positivos'

Es decir : tenemos un total de positivos (los True y los False) de :

'{:,.0f} Positivos'.format(poblacionAE+poblacionAe)
'690,525 Positivos'

entre los que son True :

'True positivos {:,.0f} de un total de {:,.0f}, representando un {:2.1f}%'.format(
    poblacionAE,poblacionAE+poblacionAe,
100*poblacionAE/(poblacionAE+poblacionAe))
'True positivos 227,850 de un total de 690,525, representando un 33.0%'

que coincide con el resultado calculado anteriormente