In [64]:
import numpy as np
import matplotlib.pyplot as pl
import openpyxl
import dateutil
import datetime
import pickle
import pystan
import seaborn as sns
import os
from hashlib import md5
import matplotlib
%matplotlib notebook
In [65]:
matplotlib.get_backend()
print(pystan.__version__)
Vamos a plantearnos como ejercicio analizar los sondeos electorales. Para ello cogemos como datos los sondeos recopilados en la página web de la Wikipedia para las elecciones de 2015 y para las de 2016. Consideramos sólo los partidos PP, PSOE, IU, Podemos y Ciudadanos. El resto de partidos se consideran metidos dentro de un mismo saco, que vienen siendo alrededor del 15% en cada sondeo.
Para modelar la evolución temporal de los sondeos para cada partido, usamos un modelado con variables latentes. La idea es que los resultados de los sondeos son sólo una estimación ruidosa y sesgada de la intención de voto real de los votantes. Por tanto, considamos que $\theta_i(t_j)$ es la subyacente intención de voto para el partido $i$ en el tiempo $t_j$. Esta variable se observa de forma ruidosa para cada empresa $k$ que realiza el sondeo, y además incluye un cierto sesgo. Por tanto, nuestro modelo de observación es:
$$ y_{ik}(t_j) \sim N(\theta_i(t_j) + b_{ik}, \sigma_{k}(t_j)^2) $$donde $b_{ik}$ es el bias que introduce la empresa de sondeos $k$ sobre el partido $i$ y $\sigma_{ik}^2$ es la varianza del sondeo realizado en el tiempo $t_j$ por la empresa $k$. Para los resultados de las elecciones, elegimos una varianza muy pequeña y no incluimos los sesgos, porque es el único momento en el que tenemos acceso directo a las intenciones subyacentes:
$$ y_{ik}(\mathrm{elecciones}) \sim N(\theta_i(t_j), \sigma_{k}(\mathrm{elecciones})^2). $$Es obvio que la suma de porcentajes de cada partido debe sumar 1 en cada sondeo, por lo cual debemos asegurarnos de que
$$ \sum_i \theta_i(t_j) = 1. $$Además, es lógico asumir que, dadas suficientes casas de sondeo, los biases sumen 0 para cada partido:
$$ \sum_k b_{ik} = 0. $$Por último, asumimos un proceso autoregresivo Markoviano para modelar la variación temporal de la intención de voto real para cada partido. Es decir, en cada sondeo, la intención de voto estará dada por una variable Gaussiana con media igual a la intención de voto anterior, y una varianza que debemos estimar para cada partido:
$$ \theta_i(t_j) \sim N(\theta_i(t_{j-1}), \sigma_i^2) $$Como es un parámetro de escala, definimos un prior poco informativo sobre $\sigma_i$:
$$ \sigma_i \sim IG(0.01,0.01). $$Por último, debemos definir un prior para $\theta_i(t_0)$, para el cual escogemos una distribución de Dirichlet con $\alpha=1$:
$$ \theta_i(t_0) \sim \mathrm{Dir}(\alpha). $$Definimos algunas funciones útiles para el parsing de la base de datos.
In [66]:
def stan_cache(model_code, model_name=None, **kwargs):
"""Use just as you would `stan`"""
code_hash = md5(model_code.encode('ascii')).hexdigest()
if model_name is None:
cache_fn = 'cached-model-{}.pkl'.format(code_hash)
else:
cache_fn = 'cached-{}-{}.pkl'.format(model_name, code_hash)
try:
sm = pickle.load(open(cache_fn, 'rb'))
except:
sm = pystan.StanModel(model_code=model_code)
with open(cache_fn, 'wb') as f:
pickle.dump(sm, f)
else:
print("Using cached StanModel")
return sm.sampling(**kwargs)
def toenglish(s):
spanish = ['ene', 'abr', 'ago', 'dic', 'de mayo de']
english = ['jan', 'apr', 'aug', 'dec', 'may']
for (j, month) in enumerate(spanish):
s = s.replace(month, english[j])
return s
def getPercentage(s):
if (s[0] not in ['0','1','2','3','4','5','6','7','8','9']):
return 0
else:
if (s.find('%') != -1):
return float(s.split('%')[0].replace(',','.')) / 100.0
else:
return float(s.split('\n')[0].replace(',','.')) / 100.0
def getSigma(s):
left = s.find('(')
right = s.find(')')
if (s[left+1:right] in ['?', '-']):
return 0.03
else:
return 1.0 / np.sqrt(float(s[left+1:right]))
def weeksDifference(d1, d2):
monday1 = (d1 - datetime.timedelta(days=d1.weekday()))
monday2 = (d2 - datetime.timedelta(days=d2.weekday()))
return int((monday2 - monday1).days / 7)
Leemos el fichero de datos obtenido de la Wikipedia y hacemos el parsing. Escogemos sólo aquellos sondeos en los que hay información para todos los partidos. Podríamos también escoger aquellos parciales porque el modelado Bayesiano los podría aceptar.
In [89]:
wb = openpyxl.load_workbook("data/sondeos.xlsx")
ws = wb.active
empresas = ['GAD3', 'Encuestamos', 'GESOP', 'Metroscopia', 'Celeste-Tel','Demoscopia Servicios', 'Simple Lógica', 'CIS',
'TNS Demoscopia', 'Invymark', 'NC Report', 'El Espanol', 'DYM', 'Sondaxe', 'Sigma Dos', 'Resultados de las elecciones']
empresaSondeoAll = []
sondeosAll = []
dateAll = []
sigmaAll = []
wbnew = openpyxl.Workbook()
# grab the active worksheet
wsnew = wbnew.active
wsnew.append(['Encuestador','Fecha','Muestra','PP','PSOE','PODEMOS','Cs','IU'])
otrosPartidos = ['F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'R']
nSondeos = 302
loopSondeos = 0
for i in range(130):
empresa = ws['A{0}'.format(i+2)].value
for (loop, emp) in enumerate(empresas):
if (empresa.find(emp) != -1):
empresaSondeo = loop
if (empresaSondeo == len(empresas)-1):
sigma = 0.0001
else:
sigma = getSigma(empresa)
PP = getPercentage(ws['C{0}'.format(i+2)].value)
PSOE = getPercentage(ws['D{0}'.format(i+2)].value)
IU = getPercentage(ws['E{0}'.format(i+2)].value)
PODEMOS = getPercentage(ws['P{0}'.format(i+2)].value)
CS = getPercentage(ws['Q{0}'.format(i+2)].value)
total = PP + PSOE + IU + PODEMOS + CS
otros = 1.0 - total
tmp = ws['B{0}'.format(i+2)].value
if (isinstance(tmp, datetime.date)):
date = tmp
else:
date = dateutil.parser.parse(toenglish(ws['B{0}'.format(i+2)].value.split('-')[-1].lower()))
tmp = date.year + (date.month-1.0) / 12.0
sondeo = [PP, PSOE, IU+PODEMOS, CS, otros]
if (0 not in sondeo):
loopSondeos += 1
sondeosAll.append(sondeo)
sigmaAll.append(sigma)
dateAll.append(date)
empresaSondeoAll.append(empresaSondeo+1)
print ("{0} - {1} {7} - s={8:4.2f} : PP={2:4.2f} - PSOE={3:4.2f} - UP={4:4.2f} - CS={5:4.2f} - Resto={6:4.2f}".format(i,
empresas[empresaSondeo], PP*100, PSOE*100, (IU+PODEMOS)*100, CS*100, otros*100, date, sigma*100))
wsnew.append([empresas[empresaSondeo],'{0}/{1}/{2}'.format(date.day,date.month,date.year-2000),1.0/sigma**2,PP*100,PSOE*100,(IU+PODEMOS)*100,CS*100,0.0])
wbnew.save('new.xlsx')
sondeosAll = np.array(sondeosAll)
nSondeos, nPartidos = sondeosAll.shape
nEmpresas = len(empresas)
print ("Número de sondeos incluidos : {0} de {1} posibles".format(loopSondeos,120))
# Compute week of every poll
weekAll = []
for i in range(nSondeos):
weekAll.append(weeksDifference(dateAll[nSondeos-1], dateAll[i]) + 1)
nDates = max(weekAll)
# Reverse all lists
sondeosAll = sondeosAll[::-1]
empresaSondeoAll.reverse()
weekAll.reverse()
sigmaAll.reverse()
dictionary = {'NPartidos': nPartidos, 'NSondeos': nSondeos, 'NEmpresas': nEmpresas-1, 'NDates': nDates, 'empresa': empresaSondeoAll, 'sondeos': sondeosAll,
'date': weekAll, 'sigmaSondeo': sigmaAll, 'alpha': np.ones(nPartidos)*1.0}
Leemos el modelo que hemos definido y hacemos el sampling.
In [85]:
date.year-2000
Out[85]:
In [60]:
f = open('model.stan', 'r')
model = f.read()
f.close()
In [61]:
out = stan_cache(model, model_name='elecciones', data=dictionary, chains=1)
In [62]:
thetaChain = out.extract('theta')['theta']
theta = np.percentile(thetaChain, [50.0, 50.0-68/2.0, 50.0+68/2], axis=0)
houseChain = out.extract('house')['house']
house = np.percentile(houseChain, [50.0, 50.0-68/2.0, 50.0+68/2], axis=0)
sigmaChain = out.extract('sigma')['sigma']
sigma = np.percentile(sigmaChain, [50.0, 50.0-68/2.0, 50.0+68/2], axis=0)
En esta representación vemos la variación de las estimaciones de voto en puntos para cada una de las encuestas que hemos considerado y la estimación subyacente que hemos inferido. Como se ve, la varianza de la estimación subyacente cambia con el tiempo, siendo muy pequeña en la fecha de las elecciones, marcada con una línea vertical.
In [63]:
colors = ["blue", "red", "violet", "orange", "yellow"]
labelsPartidos = ['PP', 'PSOE', 'PODEMOS+IU', 'Cs', 'Otros']
labelsEmpresas = ['GAD3', 'Encuest', 'GESOP', 'Metroscopia', 'Celeste',' Demosc. S.', 'S. Lógica', 'CIS', 'TNS Demos.',
'Invymark', 'NC Rep.', 'Netquest']
f, ax = pl.subplots(figsize=(13,8))
for i in range(5):
ax.plot(weekAll, sondeosAll[:,i], '.', color=sns.xkcd_rgb[colors[i]], linewidth=2)
ax.plot(np.arange(max(weekAll))+1, theta[0,:,i], color=sns.xkcd_rgb[colors[i]], linewidth=2)
ax.fill_between(np.arange(max(weekAll))+1, theta[1,:,i], theta[2,:,i], color=sns.xkcd_rgb[colors[i]], alpha=0.3)
ax.set_xlabel('Semana')
ax.set_ylabel('Fracción de votos')
ax.axvline(49)
pl.savefig("sondeos.png")
In [64]:
f, ax = pl.subplots()
sns.boxplot(data=thetaChain[:,-1,:], ax=ax)
ax.set_xticks(np.arange(5))
ax.set_xticklabels(labelsPartidos, rotation=90)
boxes = ax.artists
for (i, box) in enumerate(boxes):
box.set_facecolor(sns.xkcd_rgb[colors[i]])
pl.savefig("estimacionActual.png")
A continuación represento la evolución temporal del voto por cada partido político, incluyendo las incertidumbres estimadas para cada sondeo.
In [65]:
labelsPartidos = ['PP', 'PSOE', 'PODEMOS+IU', 'Cs', 'Otros']
f, ax = pl.subplots(nrows=2, ncols=3, figsize=(10,8))
ax = ax.flatten()
for i in range(5):
ax[i].plot(weekAll, sondeosAll[:,i], '.', color=sns.xkcd_rgb[colors[i]], linewidth=2)
ax[i].set_ylim([0,0.35])
ax[i].errorbar(weekAll, sondeosAll[:,i], fmt='none', yerr=sigmaAll, color=sns.xkcd_rgb[colors[i]], linewidth=2, ecolor=sns.xkcd_rgb[colors[i]])
ax[i].plot(np.arange(max(weekAll))+1, theta[0,:,i], color=sns.xkcd_rgb[colors[i]], linewidth=2)
ax[i].fill_between(np.arange(max(weekAll))+1, theta[1,:,i], theta[2,:,i], color=sns.xkcd_rgb[colors[i]], alpha=0.3)
ax[i].set_xlabel('Semana')
ax[i].set_xlabel('Porcentaje de votos')
ax[i].set_title(labelsPartidos[i])
pl.tight_layout()
pl.savefig("porPartido.png")
Represento a continuación las estimaciones de sesgo para cada una de las empresas encuestadoras, lo que muestra algún efecto curioso en Demoscopia Servicios que sería interesante analizar y entender por qué pasa. Nótese que un sesgo positivo implica que esa empresa subestima los votos y hace falta corregirlo por esa cantidad para obtener las observaciones.
In [66]:
f, ax = pl.subplots(ncols=3, nrows=2, figsize=(12,10))
ax = ax.flatten()
for i in range(5):
sns.boxplot(data=houseChain[:,:,i], ax=ax[i])
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labelsEmpresas, rotation=90)
ax[i].set_title(labelsPartidos[i])
pl.tight_layout()
pl.savefig("bias.png")
In [67]:
f, ax = pl.subplots(figsize=(10,6))
sns.boxplot(data=sigmaChain, ax=ax)
ax.set_xticks(np.arange(5))
ax.set_xticklabels(labelsPartidos, rotation=90)
boxes = ax.artists
for (i, box) in enumerate(boxes):
box.set_facecolor(sns.xkcd_rgb[colors[i]])
pl.savefig("variabilidad.png", facecolor=f.get_facecolor(), edgecolor='none')
In [80]:
f = open('elecciones2015/index.html','r')
lines = f.readlines()
f.close()
In [81]:
fecha = datetime.datetime.now()
lines[29] = """<a id="welcome-to-github-pages" class="anchor" href="#welcome-to-github-pages" aria-hidden="true"><span aria-hidden="true" class="octicon octicon-link"></span></a>Resultado {0}/{1}/{2}</h3>\n""".format(fecha.day,fecha.month,fecha.year)
theta = np.percentile(thetaChain, [50.0, 50.0-95/2.0, 50.0+95/2], axis=0)
lines[33] = """<td style='height: 50px; background: #02B3F4; color: white; font-size: 25px; border-radius: 1px 10px 0px 0px;'>{0:4.1f}-{1:4.1f}-{2:4.1f}</td>\n""".format(theta[1,-1,0]*100,theta[0,-1,0]*100,theta[2,-1,0]*100)
lines[34] = """<td style='height: 50px; background: #8E1744; color: white; font-size: 25px; border-radius: 1px 10px 0px 0px;'>{0:4.1f}-{1:4.1f}-{2:4.1f}</td>\n""".format(theta[1,-1,2]*100,theta[0,-1,2]*100,theta[2,-1,2]*100)
lines[35] = """<td style='height: 50px; background: #FF0202; color: white; font-size: 25px; border-radius: 1px 10px 0px 0px;'>{0:4.1f}-{1:4.1f}-{2:4.1f}</td>\n""".format(theta[1,-1,1]*100,theta[0,-1,1]*100,theta[2,-1,1]*100)
lines[36] = """<td style='height: 50px; background: #FF800E; color: white; font-size: 25px; border-radius: 1px 10px 0px 0px;'>{0:4.1f}-{1:4.1f}-{2:4.1f}</td>\n""".format(theta[1,-1,3]*100,theta[0,-1,3]*100,theta[2,-1,3]*100)
In [82]:
f = open('elecciones2015/index.html','w')
f.writelines(lines)
f.close()
In [83]:
os.system('./update.sh')
Out[83]:
In [ ]: