Análisis del comportamiento tuitero de tu congresista

Como ya saben, hace unos días terminó la discusión, e idas y venidas, acerca del grupo de trabajo de derechos humanos del Congreso, presidido por la congresista fujimorista Martha Chavez (@MarthaChavezK36).

La discusión degeneró tanto que llegó al tuiter. La congresista Martha Chavez anunciaba en tuiter sus planes de trabajo dentro de la comisión y respondía a uno que otro insulto tuitero. Era notable la cantidad de tuits emitidos por la congresista. Pero, fueron muchos tuits? pocos? en qué horas acostumbra tuitar la congresista?

Usando herramientas de Linux, Python y unas cuantas librerías «open source» podemos analizar el comportamiento tuitero de Martha Chavez.

Descargué del tuiter los 3200 tuits más recientes de la congresista. Para eso usé un cliente de tuiter usable desde la consola Linux.

t timeline -c -n 3200 MarthaChavezK36 > MarthaChavezK36.csv

Aquí ven parte de los tuits descargados (click para ampliar).

3200 tuits más recientes de Martha Chavez

Hice un gráfico del número de tuits por día, usando Python.

timeline de la congresista Martha Chavez

timeline de la congresista Martha Chavez

Este timeline comienza el 24 de julio. Vemos que tuvo bastante actividad el 28 de Julio, mediados de Septiembre (cuando se discutía sobre la unión civil de parejas del mismo sexo), primera y segunda semana de Octubre (en esa época se tuiteaba sobre la renuncia de Fujimori por fax), primera semana de Noviembre (cuando se armó el chongo de su elección como coordinadora del grupo de trabajo sobre derechos humanos).

Parece que su destitución del grupo de DDHH no hizo que Martha Chavez tuitee tanto como cuando se hablaba de la unión civil (muy revelador!).

Pero supongo que Martha Chavez tuitea en sus horas libres, cuando ya terminó sus horas de trabajo en el congreso, además de los fines de semana.

Podemos ver esto si usamos sus tuits para generar un «punchcard»:

python analizar_tuits.py MarthaChavezK36.csv | python punchcard.py -f punchcard_Martha_Chavez.png

horas de tuiteo de Martha Chavez

Esto es alucinante! La congresista tuitea todos los días de la semana. Tuitea a forro entre las 8 y 10 de la mañana (ni bien llega al Congreso?). Tuitea con mayor fuerza los días Viernes. El menor número de tuits a la 1:00pm hace suponer que a esa hora almuerza. Sábados y Domingos, no descansa, tuitea tanto como los días lunes. Y parece que se va a dormir a la 1:00 am. Al parecer duerme menos de 8 horas (eso no es saludable congresista!).

Este nivel de tuits emitidos por Martha Chavez es muy alto? muy bajo? Podemos hacer una comparación con un tuitero consumado, neto y nato. Comparemos con el Útero de Marita:

Este es el punchcard del utero.pe.

punchcard uterope

Vemos que, al parecer, el útero.pe tuitea menos que la congresista. Uterope tuitea muy poco los viernes, sábados y domingos (a excepción de las 9:00pm cuando tuitea con furia, debe ser que a esa hora pasan los noticieros dominicales). Qué hace el uterope los viernes y fines de semanas que no tuitea? Debe tener buena vida. También tuitea bastante los jueves.

Aqui les dejo el código necesario para hacer este tipo de análisis (?) con cualquier tuitero. Pero fíjense que el tuitero no ande borrando sus tuits ni use tuits programados ya que malograría el «análisis».

Sección geek

Código para producir el gráfico timeline y producir las fechas en formato unix, necesarias para dibujar el punchcard. El programa que hace el punchard lo saqué de aquí: https://github.com/aaronjorbin/punchcard.py


#! /usr/bin/env python
# -*- coding: utf-8 -*-
import sys
import codecs
import re
import datetime
import time
from itertools import groupby
import numpy as np
import matplotlib.pyplot as plt
import brewer2mpl
f = codecs.open(sys.argv[1].strip(), "r", "utf-8")
datos = f.readlines()
f.close()
timestamps = []
counting = []
x = []
for line in datos:
line = line.strip()
if re.search("^[0-9]{6,},", line):
line = line.split(",")
fecha = line[1]
unix_time = time.mktime(datetime.datetime.strptime(fecha, "%Y-%m-%d %H:%M:%S +%f").timetuple())
# correct for local time Lima -5 hours
unix_time -= 60*60*5
print unix_time
fecha = fecha.split(" ")[0]
my_time = datetime.datetime.strptime(fecha, "%Y-%m-%d")
if my_time not in timestamps:
timestamps.append(my_time)
counting.append(fecha)
if fecha not in x:
x.append(fecha)
# de reversa
timestamp = timestamps[::-1]
y_axis = [len(list(group)) for key, group in groupby(counting)]
# queremos color
set2 = brewer2mpl.get_map('Set2', 'qualitative', 8).mpl_colors
color = set2[0]
fig, ax = plt.subplots(1)
plt.plot(timestamps, y_axis, color=color)
plt.xticks(rotation="45")
plt.ylabel(u"Número de tuits por día")
plt.title(u'Actividad tuitera de Martha Chavez: timeline')
plt.tight_layout()
plt.savefig("timeline" + sys.argv[1].strip() + ".png")
sys.exit()

Inseguridad ciudadana para histéricos

Ministros y ex-primer ministros de este gobierno coinciden en pensar que la inseguridad ciudadana es una percepción, ilusión de la gente, producto de estados mentales histéricos y que no hay que quejarse tanto.

Lo cierto es que según los datos del INEI. El número total de delitos a nivel nacional está aumentando:

Número total de delitos. Fuente INEI http://www.inei.gob.pe/media/MenuRecursivo/Cap08005.xls

Figura 1. Número total de delitos. Fuente INEI http://www.inei.gob.pe/media/MenuRecursivo/Cap08005.xls

Estadísticas frecuentistas

Una regresión lineal nos confirma la tendencia:

Regresión lineal de número total de delitos versus año.

Figura 2. Regresión lineal de número total de delitos versus año.

La regresión lineal (y el gráfico) nos dice que conforme pasan los años ha aumentado la delincuencia (R2 = 0.67) de manera significativa (p-value = 0.008).

Se observa que entre los años 2008 y 2011 ocurrió un punto de quiebre y la delicuencia aumentó, pero no podemos apuntar con precisión en qué año comenzó esta racha de mayor número de delitos. Además que los datos no se ajustan muy bien a la línea de tendencia. Pareciera que el aumento de delitos no es lineal, parece ser exponencial! Esta incertidumbre es parte de las limitaciones de las estadísticas que estoy usando, esta corriente llamada estadísticas frecuentistas.

Estadísticas bayesianas

Pero afortunadamente existen las estadísticas bayesianas que nos pueden dar algo más de información respecto a este tema.

Estas estadísticas nos pueden ayudar a estimar en qué año aumentó la delincuencia. Por ejemplo, el promedio de delitos anuales antes de este incremento puede ser considerado como la variable \lambda_1, el promedio de delitos después del punto de quiebre pueder ser la variable \lambda_2, y el año en que ocurrió el punto de quiebre puede ser la variable tau.

Podemos estimar el rango de valores más probables que puede tener cada variable si es que usamos simulaciones de números aleatorios.

[Paréntesis]
Los estadísticos bayesianos conocen la probabilidad más alta de estos valores como probabilidad posterior. Por ejemplo, cuando tú amig@ lector(a) recibes un email (digamos de Gmail), la empresa Google tiene un software que aplica estadísticas bayesianas al contenido del mensaje. Lo que hace es buscar palabras clave que indiquen que el email recibido es spam. La probabilidad inicial que un mensaje sea Spam puede ser 0.5 (osea 50%), pero si el contenido tiene las palabras «viagra», «penis», «enlargement». Existirá una mayor probabilidad que este correo es spam, (a esta probabilidad se le llama probabilidad posterior ya que se obtiene luego de examinar la evidencia), y el software de Google lo enviará directamente a la carpeta Junk. Por eso las estadísticas bayesianas son importantes, y además las usas a diario sin darte cuenta.
[/Paréntesis]

Volviendo a nuestro problema, necesitamos ver cuáles son las probabilidades posteriores de nuestros datos de número de delitos a nivel nacional. Felizmente, el lenguaje de programación Python tiene una librería muy chévere para hacer estadísticas bayesianas. Es el paquete pymc. Entonces solo es cuestion de simular muchas veces los valores de números totales de delito antes y después del incremento, y el año de incremento en la tasa delincuencial (osea \lambda_1, \lambda_2, y tau).

Realicé una simulación de 50 mil generaciones usando una cadena Markov Monte Carlo, descarté las primeras 10 mil generaciones y dibujé los resultados:

Vemos que hay una diferencia notable entre el total esperado de delitos antes (\lambda_1) y después (\lambda_2) del incremento (casi 160 mil delitos antes y 230 mil delitos luego del incremento de la delincuencia).

También vemos que es más probable que en el 6to año (osea año 2011) ocurrió la aceleración de la delincuencia en el Perú.

Podemos combinar estas tres variables en un solo gráfico:

Este gráfico muestra los valores esperados de delito antes, después y durante la aceleración en el nivel de delicuencia. Vemos que esto ocurrió del año 2010 al 2011.

Ahora, pregunto qué acontecimiento ocurrió entre 2010 y 2011 que fue el causante del aumento de la delicuencia en nuestro país? El que sepa «que levante la mano».

Sección geek

Aquí los datos que he usado:


152516
153055
144205
151560
160848
181866
206610
254405

view raw

datos.txt

hosted with ❤ by GitHub

Auí está el código para hacer la regresión lineal en R:


library(ggplot2)
y <- read.csv("datos.txt", header=FALSE)
y <- as.vector(y[,1])
x <- 2005:2012
int <- lsfit(x,y)$coefficients[1]
slope <- lsfit(x,y)$coefficients[2]
p <- ggplot(,aes(x,y))
p + geom_point() + geom_abline(intercept=int, slope=slope) +
labs(title = "Perú: Número total de delitos por año")
summary(lm(x ~ y))
# R2 = 0.67
# p-value = 0.008

view raw

lm.R

hosted with ❤ by GitHub

Aquí el código para hacer el análisis bayesiano (usa Python, pymc, y prettyplotlib):


# -*- coding: utf-8 -*-
import prettyplotlib as ppl
from prettyplotlib import plt
import sys
import pymc as pm
import numpy as np
datos = np.loadtxt("datos.txt")
alpha = 1.0/datos.mean()
print alpha
print "alpha %f" % alpha
print "datos.mean %f" % datos.mean()
n_datos = len(datos)
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
print lambda_1.random()
print lambda_2.random()
tau = pm.DiscreteUniform("tau", lower=0, upper=n_datos)
print tau.random()
@pm.deterministic
def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
out = np.zeros(n_datos)
out[:tau] = lambda_1
out[tau:] = lambda_2
return out
observation = pm.Poisson("obs", lambda_, value=datos, observed=True)
model = pm.Model([observation, lambda_1, lambda_2, tau])
mcmc = pm.MCMC(model)
mcmc.sample(50000, 10000, 1)
lambda_1_samples = mcmc.trace('lambda_1')[:]
lambda_2_samples = mcmc.trace('lambda_2')[:]
tau_samples = mcmc.trace('tau')[:]
fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1)
plt.rc('font', **{'family': 'DejaVu Sans'})
plt.subplot(311)
plt.title(u'''Distribución posterior de las variables
$\lambda_1,\;\lambda_2,\;tau$''')
plt.hist(lambda_1_samples, histtype="stepfilled", bins=30, alpha=0.85,
normed=True)
plt.xlim([150000,250000])
plt.xlabel("valor de $\lambda_1$")
plt.subplot(312)
#ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype="stepfilled", bins=30, alpha=0.85,
normed=True)
plt.xlim([150000,250000])
plt.xlabel("valor de $\lambda_2$")
plt.tick_params(axis="both", which="mayor", labelsize=4)
plt.subplot(313)
w = 1.0/tau_samples.shape[0]*np.ones_like(tau_samples)
plt.hist(tau_samples, bins=n_datos, alpha=1, weights=w, rwidth=2.0)
plt.xticks(np.arange(n_datos))
plt.ylim([0, 1.5])
plt.xlim([0, 8])
plt.xlabel("valor de $tau$")
fig.set_size_inches(7,6)
fig.tight_layout()
fig.savefig("plot1.png")
fig, ax = plt.subplots(nrows=1, ncols=1)
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_datos)
for day in range(0, n_datos):
ix = day < tau_samples
expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
+ lambda_2_samples[~ix].sum()) / N
anhos = ["2005","2006","2007","2008","2009","2010","2011","2012"]
plt.plot(range(n_datos), expected_texts_per_day, lw=4, color="#E24A33",
label="expected number of text-messages received")
plt.xlim(0, n_datos)
plt.xticks(np.arange(n_datos) + 0.4, anhos)
plt.xlabel(u'Años')
plt.ylabel(u'Número esperado de delitos')
plt.title(u'''Cambio en el número esperado de delitos por año''')
plt.ylim(0, 300000)
plt.bar(np.arange(len(datos)), datos, color="#348ABD", alpha=0.65)
#plt.legend(loc="upper left")
fig.savefig("plot2.png")