Candidatos y sigmas

Vamos a hacer lo de ayer, pero introduciendo un cambio.

En vez de mirar las graficas, miramos el ajuste y sacamos una lista de candidatos que se desvian «más de un sigma» (con un sigma normalizado a ojimetro). En realidad se va a ver lo mismo que en los histogramas, pero asi me lo apunto para recordar como se hace numericamente. Fragmento de codigo:

  # a partir de aqui el analisis
  n=80
  for x in voto:
    resultado[x]+=n
    resultadoBis[x]+=1
    n=n-1
plancha={}
import numpy as np
from math import sqrt
#from itertools import accumulate #this is np.cumsum
from scipy import stats
from scipy.optimize import curve_fit
def weib(x, *p):
    N, c, s = p
    return N*stats.frechet_r.pdf(x,c,scale=s,loc=0) #N* c * y**(c-1) * np.exp(-y**c)
def gamma(x,*p):
    N, a, c, s = p
    return N*stats.gengamma.pdf(x,a,c,scale=s,loc=0) #abs(c) * x**(c*a-1) * exp(-x**c) / gamma(a)

for e in rangos:
  plancha[e]=[resultadoBis[x] for x in rangos[e]]
#  for n,x in enumerate(rangos[e])): 
#     print n,x, resultadoBis[x],nombre[x],plancha[e][n]
  xdata=range(1,len(plancha[e])+1)
  ycdf=np.cumsum(plancha[e])
  print e,":"
  resW, pcov = curve_fit(weib, xdata,plancha[e],[100000,0.7,200],maxfev=10000)#max_nfev=100000,#method="dogbox",
                             #sigma=pesos,absolute_sigma=False)
                             #bounds=[[0,0,0],[+np.inf,1,+np.inf]])
  resG, pcov = curve_fit(gamma, xdata,plancha[e],[1000,1,1,60],maxfev=100000)#,method="dogbox",
  sqdevW,sqdevG=0.0,0.0
  for n,x in enumerate(rangos[e]):
     sqdevW+=(resultadoBis[x]/weib(n+1,*resW)-1)**2
     sqdevG+=(resultadoBis[x]/gamma(n+1,*resG)-1)**2
  sqdevW=sqrt(sqdevW/len(rangos[e]))
  sqdevG=sqrt(sqdevG/len(rangos[e]))
  for n,x in enumerate(rangos[e]):
     if abs((resultadoBis[x]/weib(n+1,*resW)-1)/sqdevW) > 0.9:
       print n+1,nombre[x]," %+.2f"%((resultadoBis[x]/gamma(n+1,*resG)-1)/sqdevG)

Resultado de la ejecucion:

PeM :
61 Alberto San Juan Guijarro  +5.40
62 Jose María Gonzalez Santos  +5.01
PPT :
2 Irene Montero Gil  -1.15
4 Gloria Elizo Serrano  -1.94
5 Pablo Echenique  +2.37
6 Noelia Vera Ruiz Herrera  -1.13
9 Julio Rodríguez Fernandez  +2.36
11 Juan Manuel Del Olmo  -1.09
15 Diego Cañamero Valle  +2.24
17 Alberto Rodríguez Rodríguez  +3.81
21 Javier Sánchez Serna  -0.99
23 Pablo Manuel Fernández Alarcón  -1.03
25 Guillen Carroza Armendariz  -1.03
33 Pedro Antonio Honrubia Hurtado  -1.03
34 Rita Bosaho Gori  +2.22
37 Federico Severino Vitantonio  -1.07
RlI :
3 Juan Pedro Yllanes Suarez  -1.01
5 Pablo Bustinduy Amador  +1.31
6 Tania González Peñas  +0.91
13 Miguel Vila Gómez  +1.54
61 Pepe Viyuela  +6.64
PeEquipo :
10 Francisco Javier Fernández Ágreda  -1.11
11 Maria Teresa Rios Jimenez  +0.93
14 Jorge Gonzalez Andreu  -1.26
16 Alberto Manuel Argudo Diaz  -1.76
18 Ángel García Pageo  -1.54
20 José Juan Martínez Araque  +1.06
22 Moisés Macías Díaz  -1.31
26 Soledad Gonzalez Domingo  +1.69
27 Rosario Martín Venegas  +2.97
Blanco :
3 Maria Del Rocio Aragon Lozano  -1.06
4 Ana Garcia Arberas  +2.71
14 Jose Luis Vicien Garcia  -0.98
15 Lorena Guzmán Trochez  +1.94

La lista de Pablo Iglesias parece salir algo mas movida debido precisamente a que el ajuste trata de encajarlo. Una alternativa -que se podria aplicar para los dos o tres primeros de cada lista- es excluirlo del fit y dejar que aparezca como outlier:

PPT :
1 Pablo Iglesias Turrión  +3.55
4 Gloria Elizo Serrano  -1.60
5 Pablo Echenique  +2.26
6 Noelia Vera Ruiz Herrera  -1.08
9 Julio Rodríguez Fernandez  +1.95
11 Juan Manuel Del Olmo  -1.18
15 Diego Cañamero Valle  +1.84
17 Alberto Rodríguez Rodríguez  +3.29
19 Jose Maria Guijarro  -0.90
21 Javier Sánchez Serna  -1.00
23 Pablo Manuel Fernández Alarcón  -1.02
25 Guillen Carroza Armendariz  -1.00
33 Pedro Antonio Honrubia Hurtado  -0.92
34 Rita Bosaho Gori  +2.09
37 Federico Severino Vitantonio  -0.94
52 Marlene Pérez Sánchez  +0.93
60 Noemí Santana Perera  +0.94
61 Isabel Franco Carmona  +0.95

Una nota sobre el codigo: Estas graficas se suelen encajar a Weibull, pero en realidad el unico motivo por el que se usa esta distribucion en vez del modelo gamma es porque tiene integral exacta. En modelo numerico podriamos usar la segunda. En cualquier caso para lo que estabamos, que era un ajuste a minimos cuadrados, igual nos daba, y en el programa he probado de varias maneras. Otro asunto seria si quisieramos extraer parametros en serio, entonces habria que usar algun estimador de «likeliness» y tener cuidado con el hecho de que la distribucion esta «censurada a la derecha» por el maximo de candidatos posibles. El fit de minimos cuadrados es una forma cutre de saltarse este problema.

Extraer parametros seguramente daria una medida de lo bien que se ha confeccionado la cabecera de lista (dando un pico mas o menos grande) y del entusiasmo general que una lista da al votante que se esta dando el trabajo de decidir, en cada paso, si sigue añadiendo cruces o cambia de lista o entrega ya el voto.

Para finalizar, podemos hacer una «prediccion»: ¿Cuantos habrian votado al candidato 100, 500, o 1000 de la lista?

for x in (100,500,1000):
print x,gamma(x+1,resG),weib(x,resW)

PeM :
100 5366.80989851 5398.39128153
500 2616.94382139 2632.51790635
1000 1917.27371378 1932.12101899
PPT :
100 20493.6213224 18251.0886476
500 6462.31099485 1032.07796902
1000 3457.3992571 53.33816224
RlI :100
15909.6253579 16039.4525369
500 9013.25755656 9182.81889563
1000 7031.14078682 7221.6216423

vemos como el modelo weibull y el gamma siguen coincidiendo para 100, pero en el caso de PPT comienzan a discrepar para 500 y 1000.