Estadística

Tuve que saltarme una diapositiva en el DataBeers de Madrid del pasado jueves.

(A propósito, aquí están las 1+20 diapositivas.)

La decimonona, de la que trata la entrada, viene a hablar de lo siguiente. Tenemos una base de datos con sujetos (ids) que hacen cosas en determinados momentos. No es inhabitual calcular la frecuencia de esos sujetos así:

select id, count(*) as freq
from mytabla
where fecha between current_date - 7 and current_date
group by id
;

Esa variable se utiliza frecuentemente ya sea como descriptor de los sujetos o como alimento de otros modelos.

Si lees estas líneas y vives en Madrid o sus alrededores, casi seguro que te interesa el DataBeers de esta tarde. Que a sus tradicionales bondades añade en esta ocasión la de que participaré y hablaré durante seis minutos y veinte segundos sobre mascotas y rebaños de modelos.

¿Nos veremos?

Una de las mayores contrariedades de estar sentado cerca de alguien que es más matemático que un servidor (de Vds., no de silicio) es que oye siempre preguntar por qué. Una letanía de preguntas me condujo a leer papelotes que ahora resumo.

Primero, unos datos:

set.seed(1234)

n <- 100

x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)

y <- 0.3 + 0.2 * x1 + 0.5 * (x2 > 0) + 0.2 * rnorm(n)

Luego, un modelo:

modelo <- lm(y ~ x1)
summary(modelo)

# Call:
#   lm(formula = y ~ x1)
#
# Residuals:
#   Min      1Q  Median      3Q     Max
# -0.9403 -0.2621  0.0420  0.2299  0.6877
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  0.55632    0.03364  16.538  < 2e-16 ***
#   x1           0.21876    0.03325   6.579 2.34e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3323 on 98 degrees of freedom
# Multiple R-squared:  0.3063,  Adjusted R-squared:  0.2992
# F-statistic: 43.28 on 1 and 98 DF,  p-value: 2.341e-09

Pocos que no entiendan cómo se han generado los datos advertirían lo malo de su especificación: hemos omitido una variable explicativa cuyo efecto ha ido a incrementar el error de manera que los tests habituales de bondad de ajuste no advierten.

Dos motivos me han tenido missing estas últimas semanas. Uno es una estancia en la Universidad de Santa Catalina del Burgo de Osma. Oportunamente ubicada en las estribaciones de la muy generosa en caldos de calidad Ribera del Duero, ha sido reconvertida a la sazón en un hotel propicio para la evasión y la agrafía.

El segundo es que en horas intempestivas he estado purgando de missings unas matrices enormes y de la, se conoce, mayor trascendencia. Es un asunto delicado, jamás bien resuelto, para el que el paquete [VIM](http://cran.r-project.org/web/packages/VIM/index.html) puede proporcionar ayuda. Sobre todo en los aspectos gráficos.

Muchos cuidamos de nuestro ordenador casi como una mascota: le ponemos un nombre (a menudo escribo desde tiramisu), le hacemos algo de mantenimiento, etc. Hay quienes, incluso, decoran sus máquinas con pegatinas.

Pero llega un momento en que hay que comenzar a tratar a las máquinas no tanto como mascotas sino como rebaños. Desde una pantalla aneja a esta en la que escribo estoy manejando un clúster de más de 200 GB y 50 núcleos distribuido en varias máquinas que ni sé dónde están. Además, solo espero que crezca. Ya no cuido de una mascota; cuido de un rebaño.

Lo que escribí hace un par de días sobre procesos puntuales, ahora me doy cuenta, podía haberse resuelto con nuestro viejo amigo glm.

Ejecuto el código del otro día y obtengo (para un caso nuevo)

          mu       alfa verosimilitud delta
    1  0.4493158 0.50000000      340.6141     1
    2  0.2675349 0.40457418      307.3939     2
    3  0.1894562 0.28917407      293.4696     3
    4  0.1495654 0.22237707      287.0784     4
    5  0.1243791 0.18079703      281.3900     5
    6  0.1142837 0.14913172      284.9227     6
    7  0.1217504 0.12150745      288.5448     7
    8  0.1214365 0.10424818      289.3282     8
    9  0.1204605 0.09148817      290.9081     9
    10 0.1315896 0.07857330      295.3935    10</code>

que significa que el parámetro óptimo es delta = 5, mu = 0.124 y alfa = 0.18.

Ahora hago

    cuantos.previos <- function(i, muestra, delta){
      indices <- Filter(function(x) x < i & x > i - delta, 1:n)
      cuantos <- sum(muestra[indices])
    }

    fit.glm <- function(delta){
      prev <- sapply(1:length(muestra),
                     cuantos.previos, muestra, delta)
      dat  <- data.frame(muestra = muestra, prev = prev)

      res.glm <- glm(muestra ~ prev, data = dat,
                     family = poisson(link = "identity"))
      c(delta, res.glm$coefficients, summary(res.glm)$aic)
    }

    res.glm <- sapply(1:10, fit.glm)
    res.glm <- as.data.frame(t(res.glm))
    colnames(res.glm) <- c("delta", "mu", "alfa", "aic")

y obtengo

Tengo una serie de datos que se parecen a lo que cierta gente llama procesos puntuales y que se parecen a los que se introducen (muuuuy prolijamente) aquí. Gráficamente, tienen este aspecto:

Sobre un determinado periodo de tiempo (eje horizontal) suceden eventos y los cuento por fecha. Pero no suceden independientemente (como si generados por un proceso de Poisson) sino que tienden a agruparse: el que suceda un evento tiende a incrementar la probabilidad de que suceda otro poco después. El proceso, en una mala traducción, se autoexcita.

Partamos el tiempo en, p.e., días y contemos una serie de eventos que suceden en ellos. Es posible que esos recuentos se distribuyan según un proceso de Poisson de parámetro $latex \lambda$, que es un valor que regula la intensidad.

Si los días son homogéneos, i.e., no hay variaciones de intensidad diaria, estimar $latex \lambda$ (por máxima verosimilitud), es tan fácil como calcular la media de los sucesos por día. Pero puede suceder que la intensidad varíe en el tiempo (p.e., se reduzca los fines de semana). O que fluctúe de cualquier manera. O que haya periodos de gran intensidad y otros de calma. Es decir, que el proceso no sea homogéneo y que $latex \lambda$ varíe en el tiempo.

[…] con algo llamado aritmética política en el siglo XVII […]

El resto (son apenas ocho hojas), tenéis que leerlo (porque os intriga, ¿verdad?) aquí.

Una red bayesiana es algo de lo que ya hablé (y que me está volviendo a interesar mucho últimamente). En esencia, es un modelo probabilístico construido sobre un grafo dirigido acíclico.

Que, a su vez, es algo parecido a

que es un grafo (obviamente), dirigido (tiene flechas) y acíclico porque siguiéndolas no se llega nunca al punto de partida. Se puede construir modelos probabilísticos sobre ellos. Basta con definir para cada nodo $latex x$ la probabilidad condicional $latex P(x|A(x))$, donde $latex A(x)$ son sus padres directos. Con estas probabilidades condicionales (y un poco de esfuerzo) se puede construir la función de probabilidad completa, $latex P(x_1, \dots, x_n)$.