El lunes que viene, día 27 de mayo, impartiré un taller de… bueno, según el programa, de Principios básicos de estadística. En realidad quiero hablar principalmente de cómo evitar incurrir en el poco conocido error de tipo III —dar la respuesta correcta al problema equivocado— y, en particular, de tres de los fenómenos que nos conducen hacia él:

• La inextricable multidimensionalidad de la realidad.
• Nuestra atávica aversión a la incertidumbre.
• El poco temor de Dios con que tomamos el nombre de la causalidad en vano.

La charla formará parte del taller con el que arrancará la competición de periodismo de datos InnovaData , coorganizada por BBVA y la Fundación Ciudadana Civio, que han tenido la gentileza de invitarme.

Así que os invito a registraros en los talleres (existe la opción de seguirlos presencialmente o en remoto) y a los más intrépidos, a competir y optar a premios de hasta 3000 euros.

Bien escondidita en las encuestas que se publican, puede encontrarse a veces una ficha técnica. Y esta suele contener una frase de esta guisa: Partiendo de los criterios del muestreo aleatorio simple, para un nivel de confianza del 95 % (que es el habitualmente adoptado) y en la hipótesis más desfavorable de máxima indeterminación (p=q=0.5), el margen de error de los datos referidos al total de la muestra es de 3.2 puntos.

(Nota: la frase está extraída de aquí y hace referencia a una encuesta en la que no se usa el muestreo aleatorio simple y en cuyo original, el nivel de confianza habitualmente aceptado —errata, supongo— figura como del 99.5 % en lugar del 95 %; pero, total, tanto da: nadie lee y a nadie le importa la liturgia de la ficha técnica).

En fin, sigamos.

La cuestión es que el otro día una colega me preguntó lo siguiente: efectivamente, el margen de error en sus datos era, creo recordar, del 4 % bajo la consabida hipótesis conservadora de p=q=0.5. No obstante, su p estimada era de 0.14, muy lejos del 0.5 de la máxima indeterminación. Y su pregunta era: ¿cuál sería el margen razonable de error en ese caso?

El resto de los datos, a continuación:

N <- 546   # número de sujetos en la población
n <- 182   # número de sujetos muestreados

p <- 0.14  # proporción estimada de sujetos X

x <- ceiling(n * p)  # número obtenido de sujetos X en la muestra

Llamemos a la variable (desconocida, aleatoria) que indica el número de individuos X en la población. Lo que hemos obtenido en la muestra es una visión indirecta de , típicamente representada como

Esta expresión muestra cómo la proporción (conocida, porque se mide sobre la muestra) de sujetos X depende de la variable de interés . Y nos interesa conocer el rango de valores de compatible con el valor observado, .

Para ello usaremos el teorema de Bayes,

y dada nuestra ignorancia a priori sobre $P(\theta)$, bien podemos suponerla uniforme (esto es, independiente de ), con lo cual

Y ahora

n.reales <- 0:N
probs <- sapply(n.reales, function(y) dhyper(x, y, N-y, n, log = FALSE))
probs <- probs / sum(probs)

calcula probs, el vector de probabilidades correspondiente a la distribución a posteriori de sobre 0:N, que tiene la siguiente pinta:

Haciendo

tmp <- n.reales[order(probs)]
tmp <- tmp[cumsum(sort(probs)) > 0.05 ]
range(tmp) / N * 100
# 10.43956 18.68132

llegamos a la conclusión de que (con nuestras hipótesis) tenemos garantizado al 95% que el valor verdadero de la proporción estará confinado en el intervalo [0.104, 0.187]: el margen de error es del 4%.

Y cambiando el valor de p, también que el margen de error en el caso de máxima indeterminación no bajaría (mucho) del 6 % (y no del 4 % que habían anunciado: fíate de las colegas, fíate).

El Banco Central Europeo publicó un estudio sobre la riqueza de los hogares europeos en abril de 2013. A partir de él, el Bundesbank publicó otro informe que subrayaba las diferencias en riqueza entre los hogares alemanes y, supongo que entre otros, los españoles.

El informe de BCE recogía la media y la mediana del patrimonio de los hogares por países (junto con otras variables adicionales, como la renta, el nivel de endeudamiento, etc.). Obviamente, las medias son superiores a las medianas en prácticamente todas esas variables. El Bundesbank, en su informe, omitía las medias y presentaba únicamente las medianas, magnitudes que contribuían a subrayar una presunta pobreza relativa de los hogares alemanes respecto a los españoles.

Los medios españoles se hicieron eco de la noticia y publicaron artículos como:

• Los españoles, ¿más ricos que los alemanes? en El Mundo, en el que se confunde media y mediana.
• Un estudio revela que los hogares españoles tienen más riqueza que los alemanes en La Voz de Galicia.
• Un estudio del Bundesbank señala que los hogares españoles son más ricos que los alemanes en la Ser
• Prestigiosos medios alemanes afirman enfadados que los españoles son más ricos que los alemanes en La Información.

y otros que se pueden encontrar por ahí. Algunos de ellos —salvo el de El Mundo, que confunde los conceptos— discuten la relación entre media y mediana y la elección torticera del segundo estadístico por parte del Bundesbank.

Ya he hablado previamente sobre medias y medianas en estas páginas. Y alguna vez, imagino, habré dejado escrito que si me pusiesen una pistola en la sien y me obligasen a resumir una distribución (p.e., de riqueza de los hogares) en un único número, utilizaría la mediana.

Pero eso solo si me pusiesen una pistola en la sien. Si pudiera razonar con quien me hace el encargo, trataría de convencerlo para que mostrase gráficamente la distribución completa mediante un histograma o técnica similar. ¿Por qué restringirse a un único número cuando con poco más esfuerzo se pueden tener todos? Tenemos una obsesión enfermiza por las cifras redondas, por olvidar que la realidad danza aleatoriamente alrededor de esas estacas cognitivas —las que conforman el temario de esos cursos de estadística introductoria— que clavamos en nuestros cerebros para autosugestionarnos con que comprendemos los fenómenos.

No nos restrinjamos innecesariamente, mostremos las distribuciones completas y, entre otras cosas, nos ahorraremos discusiones bizantinas sobre cuántos ángeles pueden danzar sobre la punta de un alfiler.

Coda: Guillermo de la Dehesa ha publicado recientemente un artículo, ¿Son los hogares españoles más ricos que los alemanes?, que aunque se atora momentáneamente en la inane discusión sobre la media y la mediana, tiene a bien ahondar en otros aspectos de las estadísticas, como lo que miden y lo que no miden realmente y que explican bastante bien las causas de la aparente paradoja que ha motivado tanta ida y venida de artículos.

El jueves 16 de mayo hablaré en el Grupo de Interés Local de Madrid de R sobre lematizadores probabilísticos.

Hablaré sobre el proceso de lematizacion y trataré de mostrar su importancia dentro del mundo del llamado procesamiento del lenguaje natural (NLP). La lematización es un proceso humilde dentro del NLP del que apenas nadie habla: su ejercicio solo ha hecho famoso a Martin Porter. Lo eclipsan otras aplicaciones más vistosas, como el siempre sobrevalorado análisis del sentimiento. Sin embargo, es una pieza fundamental que subyace (o debería subyacer) en cualquier aplicación seria que analice textos.

En la charla repasaré las tres grandes familias de soluciones para el problema de la lematización:

• las basadas en reglas duras,
• las basadas en diccionarios y, finalmente,
• las más interesantes, las probabilísticas.

Y, en particular, describiré con cierto detalle —aunque tratando de obviar los aspectos técnicos más áridos— un algoritmo que combina oportunísticamente diccionarios y modelos ocultos de Markov y que debería ver la luz en producción dentro del conjunto de APIs lingüísticas de Molino de Ideas.

Esta entrada tiene como prerrequisito las dos que la preceden: esta y esta.

Si es una muestra de una distribución de probabilidad regular y extendida, entonces es una muestra de , que es otra distribución de probabilidad

• regular (porque el logaritmo es una función creciente) y
• extendida (aunque hay que convenir que menos: el logaritmo achica los números grandes).

Por lo tanto, cabe esperar que también la parte decimal de tenga una distribución uniforme sobre el intervalo [0,1). Luego cumple la Ley de Benford (véase la condición suficiente). Esto se debe a esa (¿contraintuitiva?) propiedad del logaritmo decimal: convertir el dígito más significativo de un número, el primero, en la parte menos significativa de su logaritmo, la que sigue a la coma.

Tres notas de rigor:

• En lugar de podrían usarse otras funciones (el cuadrado, la raíz cuadrada, etc.) que también transforman distribuciones regulares y extendidas en otras que lo son igualmente. Pero se perdería la magia de la relación entre la parte fraccionaria con el primer dígito.
• La parte fraccionaria de una distribución regular y extendida es aproximadamente uniforme. La uniformidad solo se garantiza en el límite (conforme la distribución se hace más y más extendida sobre la recta real). Es posible (cuestión que exploré aquí) que los primeros dígitos de muestras de determinadas distribuciones no sigan la Ley de Benford.
• Queda ver cuáles son las razones (¿sicológicas?) que llevarían a los humanos a inventar secuencias de números que no obedecen una ley extendida y regular. En particular, que violan la regularidad.

En la última entrada de la serie abundaré esa tercera nota y hablaré de posibles extensiones que no son sino ocurrencias mías.

Sigo con mi lacónica serie sobre data.table.

La protagonista:

frases[sample(1:nrow(frases), 3),]
#pos.es pos.en length.es length.en en        es frase          tfe      qjilm          num
#1:     15     43        72        72  i        de  2632 4.881416e-02 0.01369863 6.686871e-04
#2:     33     48        46        48  X    países  5321 2.726146e-06 0.02040816 5.563563e-08
#3:      2     35        53        66 in preguntar  4582 2.424379e-08 0.01492537 3.618476e-10
dim(frases)
#[1] 6340091      10

El tiempo:

system.time({
  setkey(frases, "frase", "es")
  denominadores <- frases[, sum(num), by = key(frases)]
  setnames(denominadores, c("frase", "es", "den") )
  frases <- merge(frases, denominadores)
  frases$delta <- frases$num / frases$den
})
#user  system elapsed 
#5.628   0.208   5.841 

En particular,

system.time( denominadores <- frases[, sum(num), by = key(frases)] )
#user  system elapsed 
#0.228   0.000   0.228 

Increíble, ¿no?

El otro día tropecé con un problema de rendimiento con R y al utilizar Rprof() encontré muchas llamadas a funciones que yo no hacía directamente.

La principal sospechosa era la función daply (del paquete plyr) que parecía depender de bastantes otras. Uno puede navegar el código de las funciones para identificar esas dependencias, pero, mirad qué maravilla:

library(mvbutils)
library(plyr)
foodweb(find.funs("package:plyr"), prune = "laply")

genera

Ahí se ve la dependencia de daply con respecto a laply. Y uno adquiere, además, una visión panorámica del paquete plyr.

La función foodweb tiene argumentos adicionales (prune es uno de los más útiles) que los interesados podrán encontrar en ?foodweb.

Motivado por los experimentos de Gregorio Serrano con shiny e ilustrado por la charla que dio en el Grupo de Usuarios de R de Madrid, decidí colgar el otro día un entretenimiento que ocupó la mañana de un domingo —y las mañanas de mis domingos son proverbialmente breves— en la red.

Se trata de una aplicación que distingue el idioma en que está escrito un texto dentro de una selección de ellos: español, italiano, latín, francés, portugués y catalán.

¡Invitados quedáis a probarlo!

(Nota: esta aplicación me fue sugerida, en forma de reto, por Elena Álvarez de Molino de Ideas).

(Coda: cualquier día, por chinchar, hago uno que distinga entre catalán, valenciano y mallorquín.)

Continuamos hoy nuestra serie sobre la llamada ley de Benford discutiendo la distribución de la parte fraccionaria de las muestras de una distribución.

La parte fraccionaria de un número es, para entendernos, lo que va detrás de la coma. Técnicamente, x - floor(x). ¿Le sorprendería a alguien la parte fraccionaria de una secuencia aleatoria de números no tenga una distribución uniforme sobre [0,1)?

Obviamente, si los números son enteros no. ¿Pero si siguen la distribución normal? Se puede probar, de hecho, que si la serie sigue una distribución de probabilidad que sea

• regular, es decir, que no tenga picos extraños y, más en concreto, cuya función de densidad crezca hasta cierto punto y decrezca de él en adelante y
• extendida, es decir, que cubra un rango amplio de valores (p.e., la recta real entera),

entonces la distribución de la parte fraccionaria de sus muestras serán aproximadamente uniformes. Y lo serán tanto más cuanto menor sea el máximo de la función de distribución. La referencia, el artículo Pourquoi la loi de Benford n’est pas mystérieuse de Nicolas Gauvrit y Jean-Paul Delahaye.

Esto se verifica fácilmente en ciertos casos. Por ejemplo,

x <- rnorm(100000)
x <- x - floor(x)

par(mfrow=c(1,2))
hist(x, col = "grey")
qqplot(x, runif(100000), main = "qqplot")

que produce

En la siguiente entrega analizaremos qué tiene que ver esto con la ley de Benford.

Los protagonistas (tres tablas grandecitas):

dim(qjilm)
# [1] 3218575 5
dim(tf)
# [1] 6340091 7
dim(tfe)
#[1] 1493772 3

head(qjilm, 2)
#pos.es length.en length.es pos.en qjilm
#1 1 2 1 1 0.8890203
#2 1 2 1 2 0.1109797

head(tf, 2)
#frase es pos.es length.es en pos.en length.en
#1 996 ! 42 42 ! 43 44
#2 1231 ! 37 37 ! 37 38

head(tfe, 2)
#en es tfe
#1 ! ! 4.364360e-01
#2 ! !" 4.945229e-24

El objetivo (cruzarlas por los campos comunes):

res <- merge(merge(tf, tfe), qjilm)

El tiempo (usando merge):

res <- merge(merge(tf, tfe), qjilm)
#user system elapsed 
#442.991 2.496 446.832 

dim(res)
#[1] 6340091 9

Y con data.table:

library(data.table)

system.time({
  res.dt <- merge(data.table( tf,  key = c("en", "es")),
                  data.table( tfe, key = c("en", "es")) )

  res.dt <- merge( setkeyv(res.dt,   cols = c("pos.es", "pos.en", "length.es", "length.en")),
                   data.table(qjilm, key  = c("pos.es", "pos.en", "length.es", "length.en") )
              )
})
#user system elapsed 
#32.070 0.012 32.118 

dim(res.dt)
#[1] 6340091 9

Y, finalmente, suponiendo que los data.tables ya tienen asociado un índice de antemano:

tf.dt  <- data.table( tf,  key = c("en", "es"))
tfe.dt <- data.table( tfe, key = c("en", "es"))

system.time( res <- merge(tf.dt, tfe.dt) )
#user system elapsed 
#3.464 0.000 3.466 

dim(res)
#[1] 6340091 8

Resumen:

• Había hecho unas pruebas con data.table previamente que no resultaron del todo satisfactorias.
• Anoche, data.table me sacó de un apuro muy serio.
• Ahora soy fan.
• Gracias a data.table, el límite de tamaño de los conjuntos de datos con los que soy capaz de trabajar razonablemente con R ha crecido en todo un orden de magnitud: ya no me asusta que me hablen de millones.