Estadística

Diferencia de medias a la bayesiana con salsa de stan

El habitual problema de la diferencia de medias suele formularse de la siguiente manera: hay observaciones $latex y_{1i}$ e $latex y_{2i}$ donde

$$ y_{ji} \sim N(\mu_j, \sigma)$$

e interesa saber si $latex \mu_1 = \mu_2$. Obviamente, se desconoce $latex \sigma$. De cómo resolvió Gosset el problema están los libros de estadística llenos. En R,

set.seed(1234)
N1 <- 50
N2 <- 50
mu1 <- 1
mu2 <- -0.5
sig1 <- 1
sig2 <- 1
 
y1 <- rnorm(N1, mu1, sig1)
y2 <- rnorm(N2, mu2, sig2)
 
t.test(y1, y2)
# Welch Two Sample t-test
#
# data:  y1 and y2
# t = 4.7059, df = 95.633, p-value = 8.522e-06
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   0.5246427 1.2901923
# sample estimates:
#   mean of x  mean of y
# 0.5469470 -0.3604705

En rstan,

Ver 53000 filas

Me preguntaban cómo ver con R una tabla con 53000 filas. Mi yo menos diplomático quiso contestar: define ver. Lo reformulé más amablemente y se me contestó: como en Excel.

La pregunta es: ¿permite Excel ver 53000 registros? De hecho, ¿se pueden ver 53000 registros? Impresos a razón de línea por centímetro, ocuparían 530 metros y andar a paso vivo del primero al último costaría cinco minutos.

Con 53000 registros, ver (como trasunto de entender) es una cosa distinta de tener delante. Lo siento, pero ver otra cosa que la facturación de los últimos quince días o los movimientos de la cuenta del último mes es algo distinto de lo que vacuamente promete Excel.

Cartogramas vs huertogramas

Esto es un huertograma:

huertograma_es

Tiene la propiedad de que casi todos los pixels están encima de un huerto (o un erial, o en un cerro,…).

Este es otro huertograma:

huertograma_uk

Y esta es la misma información (resultados de las elecciones de 2015 en el RU) sobre un fabuloso cartograma:

cartograma_uk

¿Os habéis fijado cómo esa casi indistinguible mancha roja en la zona de Londres del huertograma adquiere su debida relevancia en el cartograma?

Banzhaf y las elecciones que se nos vienen

Es pertinente rescatar una entrada de hace tres años sobre D’Hondt y Banzhaf. En el enlace, los detalles.

Me limitaré a actualizar el código de la función para que muestre las alianzas (algunas enteramente esperpénticas) posibles, que queda de la forma

banzhaf <- function(x){
  x <- -sort(-x)
  x <- x/sum(x)

  foo <- function(a,b,p){
    if(p>1/2)
      return(list(a))

    if (length(b)==0)
      return(NULL)

    b.prima <- b[-1]
    delta <- b[1]
    p.delta <- x[delta]

    return(c(foo(c(a,delta), b.prima, p+p.delta), foo(a,b.prima,p)))
  }

  res <- foo( NULL, names(x), 0)
  print(res)
  sort( table(unlist(res)) / length(res) )
}

y a aplicarlo sobre algunos casos de la más rabiosa actualidad que Leda Duelo ha tenido la gentileza de preparar para mí y, a través de esta página, para ti también. Son los que siguen.

Frecuencias naturales (y consumo de cerveza)

Las frecuencias naturales se utilizan como alternativa a los porcentajes para expresar probabilidades en lugar de, por ejemplo, porcentajes.

frecuencias_naturales

El gráfico anterior está extraído de este documento en el que sus autores argumentan que transmite más eficazmente la idea de probabilidad que los porcentajes desnudos tan habituales.

Entienden que es preferible decir que de cada 100 litros de cerveza vendidos en España, 20 se distribuyen en botella, 30 en lata y 30 en barril (¡ya sé que no suman 100!) que reescribir la información anterior en forma de porcentajes. Eso, sí, respetando una misma cantidad de partida y porsupuestísimo, no escribiendo, como aquí, que

Cualquier parecido con la realidad es pura coincidencia

@adolflow (en persona) viene hoy y me dice si lo he visto. ¿Qué cosa? Se refiere a lo que han publicado en El Español, España en Cifras. Lo miro por encima y encuentro

tasa_paro_municipal

¡Tasa de paro municipal! Lo siento, @adolflow, pero tal cosa no existe. No, no es que los datos sean secretos, no sean transparentes, no sean reutilizables. Es, simplemente, que no existe.

¿Peros?

No, no hay peros. Fijáte: hay 8000 municipios y la EPA se basa en una encuesta de unos 60000 hogares. ¡Echa cuentas!

Pocos de los encuestados...

Como aragonés, a veces me interesa el estado de ese idioma que algunos quieren convencerme de que me es propio. En la Wikipedia hay un mapa que indica la presunta distribución de las distintas lenguas en Aragón y tienen marcado de rojo zonas que no conozco mal y en las que jamás he oído hablar en tal cosa.

Fuera de los mapas que se colorean ateniéndose a criterios poco transparentes, ¿qué nos dicen los estudios serios que puedan haberse hecho sobre los hablantes de esa lengua? Uno de los estudios más recientes que he visto (2006), Usos del aragonés en el Aragón aragonesparlante, en la página 95 y siguientes de esto, describe los resultados de una encuesta que realizaron sus autores a una muestra de 431 sujetos (n = 431) de 16 y más años residentes en los municipios de la zona incluida en el dominio lingüístico del aragonés.

A cuento de mi clase práctica de modelos no supervisados

A cuento de la sesión práctica de modelos no supervisados que impartiré este sábado y que estoy preparando justo ahora, traigo a la atención de mis lectores una imagen que el asunto me sugiere:

dog2

La fuente también vale la pena. Aunque habla de otra cosa.

Nota: releyendo la entrada antes de publicarla definitivamente, me doy cuenta de que igual estoy siendo excesivamente sutil.

Intervalos de credibilidad para la beta: una alternativa

A partir de los comentarios de Olivier Núñez a mi entrada anterior casi homónima, se nos ha ocurrido a ambos de forma independiente y simultánea una manera alternativa de calcular el intervalo: minimizando su longitud.

a <- 3
b <- 5
alfa <- 0.05

# versión de la entrada anterior:
f <- function(x){
  (dbeta(x[2], a, b) - dbeta(x[1], a, b))^2 +
    (pbeta(x[2], a, b) - pbeta(x[1], a, b) -1 +  alfa)^2
}

res <- optim(c(a/(a+b), a/(a+b)), f)
res$par
#[1] 0.08052535 0.68463436

# nueva versión
f.alt <- function(x){
  qbeta(x+0.95, a, b) - qbeta(x, a, b)
}

res.alt <- optim(0.025, f.alt)
qbeta(c(res.alt$par, res.alt$par + 0.95), a, b)
#[1] 0.08054388 0.68464900

Intervalos de credibilidad para la distribución beta

Tengo un parámetro, la p de una binomial, que supongo distribuido según una beta. Me da igual para el caso si la distribución a priori es o no informativa. Solo digo que la distribución a posteriori es otra beta con parámetros a y b.

Quiero construir un intervalo de credibilidad para p, es decir, encontrar un subintervalo de [0,1]

  • dentro del cual la densidad de la beta sea mayor que fuera y que
  • capture $latex 1-\alpha$ de la probabilidad total.

Gráficamente,

¿13.100 más/menos cuántos parados menos?

¿Cuál es la cifra de variación del número de parados de la que hablan la última EPA y los medios? 13100.

¿Más menos cuánto? Según el INE, el error de muestreo relativo, $latex \sqrt{V(\hat{\sigma}}$ a nivel nacional en términos porcentuales es

error_relativo

Es decir, el intervalo de confianza para la cifra de parados tendría una anchura como de 100k sujetos. Obviamente, eso impide calcular variaciones de un orden de magnitud menor.

¿Si un día faltan 21.63 euros en caja?

Si un día faltan 21.63 euros en caja se cuenta y se recuenta. Se revisan los tiques, se comprueban los pagos con tarjeta, se vuelven a sumar los pagos a proveedores, etc. Hasta que, con suerte, alguien encuentra algo y la diferencia se reduce a, digamos, 3.92 euros. Pero cuando la diferencia es de 2.15… se da por buena sin más.

Cuando el t-test da un p-valor de .058, se revisan los números, se reestudia la carga y manipulación de datos, se replantea si el caso 194 es o no un outlier, etc. Pero si el p-valor es 0.036, nada de eso ocurre. Nadie revisa caso 194. ¡Ni falta que hace!

La información es sorpresa

Hace unos días publiqué esto en Twitter:

David Cabo, muy oportunamente, denunció

Cosa que no niego. La frase que resumía el enlace tiene esa pintaza. No obstante, el artículo al que apunta es una elaboración de esa frase. El artículo, además, incluye (y no es habitual) referencias a dos artículos académicos (que no he consultado) que, entiendo, tratan y desarrollan la cuestión.