Simulación

El otro día publiqué un pequeño fragmento de código,

a <- 2.89
b <- 36.81
sample_dist <- function() rbeta(1, a, b)

sample_p <- function(y){
  candidate <- sample_dist()
  my_sample <- runif(1)

  if (y == 1) if (my_sample < candidate) return(candidate)
  if (y == 0) if (my_sample < 1 - candidate) return(candidate)

  sample_p(y)
}

p1 <- replicate(100000, sample_p(1))
p0 <- replicate(100000, sample_p(0))

auc <- mean(p1 > p0)
auc

que había usado antes aquí, para muestrear unas distribuciones relacionadas con el cálculo del AUC en modelos perfectamente calibrados. Lo había escrito meses atrás y supongo que me pasó como a la mayoría de mis lectores: darlo por bueno primero y usarlo después suponía todo un acto de fe (en mí, además).

[Menos mal que se me ha ocurrido buscar en mi propio blog sobre el asunto y descubrir —no lo recordaba— que ya había tratado el asunto previamente en entradas como esta, esta o esta.]

El problema de la propagación de errores lo cuentan muy bien Iñaki Úcar y sus coautores aquí. Por resumirlo: tienes una cantidad, $X$ conocida solo aproximadamente —en concreto, con cierto error— e interesa conocer y acotar el error de una expresión $f(X)$.

Está extraído de aquí y dice los siguiente:

Las Proyecciones de Población constituyen una simulación estadística de la población que residiría en España, sus comunidades autónomas y provincias en los próximos años, así como de la evolución de cada uno de los fenómenos demográficos básicos asociados, en caso de mantenerse las tendencias y comportamientos demográficos actualmente observados.

Para interpretar correctamente los resultados de las Proyecciones de Población es importante distinguir entre previsiones y proyecciones demográficas. Si bien pueden emplear el mismo método de cálculo, difieren en la filosofía.

Esta entrada es casi una referencia para mí. Cada vez tiro más de lme4 en mis modelos y en uno en concreto que tengo entre manos toca simular escenarios. Para lo cual, simulate.merMod.

Véamoslo en funcionamiento. Primero, datos (ANOVA-style) y el modelo que piden a gritos:

library(plyr)
library(lme4)

a <- c(0,0,0, -1, -1, 1, 1, -2, 2)
factors <- letters[1:length(a)]

datos <- ldply(1:100, function(i){
    data.frame(x = factors, y = a + rnorm(length(a)))
})
modelo <- lmer(y ~ (1 | x), data = datos)

El resumen del modelo está niquelado:

summary(modelo)

# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ (1 | x)
# Data: datos
#
# REML criterion at convergence: 2560.3
#
# Scaled residuals:
#     Min      1Q  Median      3Q     Max
# -3.6798 -0.6442 -0.0288  0.6446  3.3582
#
# Random effects:
#     Groups   Name        Variance Std.Dev.
# x        (Intercept) 1.5197   1.2328
# Residual             0.9582   0.9789
# Number of obs: 900, groups:  x, 9
#
# Fixed effects:
#     Estimate Std. Error t value
# (Intercept) -0.009334   0.412212  -0.023

En particular,

Por referencia y afán de completar dos entradas que hice hace un tiempo sobre el método delta, esta y esta, dejo constar mención al paquete propagate, que contiene métodos para la propagación de la incertidumbre.

Para desavisados: si $x \sim N(5,1)$ e $y \sim N(10,1)$, ¿cómo sería la distribución de $x/y$? Etc.

NIMBLE ha sido uno de mis más recientes y provechosos descubrimientos. Mejor que hablar de él, que otros lo harán mejor y con más criterio que yo, lo usaré para replantear el problema asociado el método delta que me ocupó el otro día.

Casi autoexplicativo:

library(nimble)

src <- nimbleCode({
    T_half <- log(.5) / k
    k ~ dnorm(-0.035, sd = 0.00195)
})

mcmc.out <- nimbleMCMC(code = src,
    constants = list(),
    data = list(), inits = list(k = -0.035),
    niter = 10000,
    monitors = c("k", "T_half"))

out <- as.data.frame(mcmc.out)

# hist(out$T_half), sd(out$T_half), etc.

Cosas:

Tenemos en Circiter un proyecto sobre el que no puedo dar muchos detalles, pero que vamos a plantear (en versión muy resumida) como un modelo que alimenta una simulación.

El modelo no va a ser un modelo sino un modelo por sujeto (rebaños, los llamamos aquí). Los modelos serán, casi seguro, modelos mixtos (lmer/glmer).

Pero claro, si usas un modelo, por muy mixto que sea, con intención de simular, predict se queda muy corto (¡siempre da la el mismo resultado!).

Debo esta entrada a la diligencia de Juanjo Gibaja, que se tomó la molestia de ubicar los teoremas relevantes en el libro Simulation and the Monte Carlo Method de Rubinstein y Kroese.

Esencialmente, como la distribución normal multivariante (con matriz de covarianzas I) es simétrica, entonces, dadas $X_1,\dots, X_m \sim N( 0, I_n )$ independientes, los m puntos del espacion n-dimensional $X_i/| X_i |$ siguen una distribución uniforme sobre su esfera (su superficie, vale la pena reiterar) unidad.

El otro día leí el artículo A Prototype Model of Stock Exchangede G. Caldarelli, M. Marsili y Y.C. Zhang. La promesa que me ofrecía era la de la creación de un sistema relativamente realista de los agentes que operan en los mercados financieros que diese lugar a una evolución de precios con propiedades similares a las observadas.

Sin embargo, el planteamiento, interesante en un principio, se deshinchó enseguida:

• El modelo planteado por los autores ni siquiera aspira a representar los aspectos más distintivos del mercado: en lugar de agentes tremendamente desiguales en tamaño y entrelazados en una maraña de dependencias e influencias mutuas, los agentes son todos equivalentes en tamaño (si bien es cierto que en el estado estacionario de la simulación los ingresos adquieren una distribución dada por una ley de potencias) y que actúan de manera independiente entre sí una vez observados los precios en el mercado.
• Los resultados, una serie temporal de precios, es calificada por los autores como muy rica, aunque enseguida pasan, en un dechado de honradez, a apuntar diferencias más o menos manifiestas entre sus características estadísticas y las observadas en mercados reales.

Entiendo y aplaudo el virtuosismo técnico empleado por los autores del artículo y la implementación de los algoritmos involucrados. No obstante, tras leerlo, me embriaga una extraña sensación que no debe de ser muy distinta de aquellos comensales a los que se les anunció paella, se les advirtió que no traía arroz ni gambas y comprobaron después que sabía a judías de lata.